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•  Three-dimensional,  non-unifonn,  axisymmetric  medium 

•  Arbitrary  particle  mixtures  and  size  distribution 

•  Line-by-line  and  band-model  treatments 

•  Alternate  levels  of  engineering  approximation 

These  requirements  are  met.  by  developing  two  basic  approaches,  a  6-flux 
and  a  N-flux  approach.  This  distinction  refers  to  the  associated  sub¬ 
division  in  angular  space,  six-orthogonal  coordinates  versus  N  non- 
orthogonal  coordinates  in  azimuth  and  elevation,  thus  providing  the  alter¬ 
natives  of  computational  simplicity  of  the  S-flux  model  with  generality  of 
the  N-flux  model.  Numerous  parametric  calculations  were  performed  to  study 
the  effect  of  optical  depth,  scattering  albedo,  phase  function  and  medium 
geometry  on  the  target  signature  and  computational  requirements  .^Comparisons 
with  independent  calculations  such  as  the  two-oimensional  Monte 'Carlo 
method  were  made.  The  results  show  the  different  degrees  of  accuracy 
achievable  with  the  various  models  as  functions  of  the  radiation  field 
parameters.  To  the  authors'  Knowledge,  the  N-flux  model  represents  the 
first  application  of  general  finite-element  techniques  to  emitting  cylin¬ 
drical  media.  It  provides  a  powerful  tool  for  the  analysis  of  plume 
scattering  problems  of  arbitrary  complexity.  The  accuracy  of  the  prediction 
is  limited  only  by  the  computational  resources  of  the  user. 

The  solution  of  the  equation  of  radiative  transfer  is  derived  on  a  mono¬ 
chromatic  basis.  For  the  application  to  the  finite  spectral  intervals 
associated  with  the  use  of  instruments  with  finite  resolution,  spectral 
averaging  procedures  are  required.  Thus, an  extension  of  earlier  band  model 
formalisms  for  regular  and  random  arrays  of  spectral  lines  -was  developed 
to  include  scattering  as  well  as  absorption  in  a  medium  in  which  radiative 
transfer  is  described  by  the  two-flux  method.  Representations  of  the 
averaged  transmissivity  for  the  Elsasser  model  and  the  random  model  with 
exponential-tailed  S" 1  line  intensity  distribution  were  derived  and 
numerical  results  were  obtained  for  large  ranges  of  optical  depths  of 
absorbing  and  scattering  media.  The  effects  of  inhomogeneities  away  from 
the  line-of-sight  are  accounted  for  by  a  contribution  to  the  source 
function  which  is  calculated  by  a  simplified  six-flux  technique. 

A  survey  of  available  optical  parameters  (complex  index  of  refraction  as 
a  function  of  wavelength  and  temperature)  was  made  for  particles  of  interest. 
The  compilation  by  Aerospace  was  found  useful  as  a  starting  base  for  a 
number  of  relevant  species,  namely  carbon,  aluminum  oxide,  magnesium  oxide 
and  zirconium  oxide.  Data  for  other  species  must  yet  be  obtained.  The 
complex  indices  of  refraction  for  carbon  and  aluminum  oxide  were  used, 
together  with  postulated  size  distributions,  to  calculate  the  volume  ex¬ 
tinction,  scattering  and  absorption  coefficients  and  the  angular  scattering 
phase  function  by  the  Mie  Theory  for  homogeneous,  spherical  particles  in 
the  l-25|jm  wavelength  region.  Based  on  these  results,  a  tentative  plan 
was  developed  for  the  organization  of  the  Data  File  and  its  integration  into 
the  overall  code  architecture. 
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INTRODUCTION 

Validated  models  to  predict  rocket  plume  IR  signatures  must 
be  capable  of  treating  a  wide  range  of  propulsion  and  trajectory 
parameters  as  well  as  handling  a  number  of  sensor-related  vari¬ 
ables.  While  significant  advances  have  beet*  made  in  the  state- 
of-knowledge  of  plume-related  radiative  transport  phenomenologies 
during  the  last  tan  years,  particularly  in  regard  to  advanced 
liquid  propellants,  prediction  methods  have  not  been  developed  to 
adequately  treat  the  coupling  of  gas  and  particle  radiative 
transfer  in  the  plumes  of  solid  propellant  missiles.  It  is 
particularly  important  that  models  be  developed  to  accurately 
treat  plumes  containing  scattering  particles  while  at  the  same 
time  also  yielding  accurate  results  for  conventional  liquid  pro¬ 
pellant  plumes. 

The  primary  goals  of  the  SIRRK  program  are  (i)  to  develop 
an  analytical  methodology  to  treat  coupled  gas/particle  radiative 
transport  effects  and  (ii)  to  integrate  this  analytical  capability 
together  with  the  state-of-the-art  treatments  in  hot  gas/ 
atmospheric  radiative  transport  (i.e.,  band  model  methodologies, 
band  model  parameters,  line-by-line  methodologies,  spectral  line 
parameters)  into  a  comprehensive,  user-oriented  plume  radiation 
model . 
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State-of-the-art  plume  signature  predictions  require  four 
basic  steps.  First,  the  species  concentrations  and  temperature 
in  the  exit  plane  are  calculated  through  a  combustion/nozzle 
code.  Second,  these  exit  plane  parameters  are  input  to  a  gas- 
dynamic  code  which  calculates  the  exhaust  plume  flow  field. 

Third,  the  flow  field  serves  as  input  to  a  radiation  code  that 
calculates  the  intrinsic  plume  signature.  Finally,  the  apparent 
plume  signature  viewed  by  a  remote  observer  sensor  is  obtained 
through  attenuating  the  source  intensity  by  the  atmospheric 
absorption.  When  the  source  radiation  is  correlated  with  the 
atmospheric  attenuation  (radiating  plume  species  are  also  present 
in  the  atmosphere  as  absorbers),  the  last  two  steps  are  combined 
into  a  single  calculation. 

The  present  study  focuses  on  the  third  and  fourth  component , 
i.e.,  the  calculation  of  the  intrinsic  and  apparent  plume  signature 
The  program  is  divided  into  four  major  tasks.  Task  I  calls  for 
the  formulation  of  the  SIRRM  code  methodology.  In  Task  II,  the 
code  is  to  be  developed  and  its  capabilities  demonstrated.  The 
validation  and  application  of  SIRRM  takes  place  in  Task  III.  The 
documentation  and  code  acceptance  is  done  in  Task  IV. 

Task  I  is  concerned  with  the  development  of  an  integrated 
physical  model  which  accounts  for  the  essential  spectroscopic 
and  geometric  effects  governing  single  and  multiple  photon 
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interactions  in  a  two-phase  gas-particle  mixture.  A  key 
element  in  this  task  is  the  theoretical  description  of  coupling 
phenomena  associated  with  the  spectral  and  spatial  redistribution 
of  radiation  by  scattering.  In  order  to  provide  a  correct 
description,  the  model  must  account  for  the  interaction  of 
emitted,  absorbed,  and  scattered  radiation  from  line  and  continuum 
sources  along  multiple,  intersecting  paths,  The  task  is  divided 
into  two  subtasks.  Subtask  IA  deals  with  the  development  of 
the  coupled,  multi-scattering  model,  while  Subtask  IB  calls  for 
the  formulation  of  the  overall  methodology,  i.e.,  the  integration 
of  the  scattering  model  with  the  other  parts  of  the  code.  This 
document  is  an  interim  technical  report  describing  the  development 
of  the  coupled  multi-scattering  model  performed  under  Subtask  IA. 

Section  2  contains  the  description  of  our  efforts  in  deve¬ 
loping  the  multi-scattering  model  with  cylindrical  geometry, 
using  both  engineering  approximations  (six-flux)  and  exact 
solutions  (N-flux) .  Diagnostic  calculations  with  these  models 
and  comparisons  of  the  results  with  those  obtained  by  independent 
codes  indicate  that  the  required  objectives  have  been  met. 

Since  the  present  scattering  model  is  developed  for  the 
monochromatic  case,  a  parallel  study  was  begun  to  treat  the 
spectral  averaging.  The  step  forward  to  a  band-model  formulation 
coupling  particle  scattering  with  spectral  line  absorption  is 
taken  in  Section  3,  where  the  Elsasser  and  random  model  formulations 
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are  developed  for  a  two-flux  scattering  model  and  conceptual 
approaches  for  the  six-flux  and  N-flux  models  are  outlined. 

The  fundamental  inputs  into  any  scattering  models  are 
the  optical  parameters  of  the  particles.  These  are  discussed 
in  Section  4,  including  the  Mie  theory,  the  indices  of  refraction 
and  absorption,  the  scattering  due  to  single  particles  and 
particle  size  distributions,  and  the  handling  of  the  Data  File 
as  it  would  be  integrated  into  the  overall  code  architecture. 
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SCATTERING  MODEL 


2 . 1  Introduction 

The  objective  of  Task  IA  is  to  formulate  a  scattering  model 
that  incorporates  the  following  capabilities: 

•  Coupled  emission,  absorption  and  scattering  treatment 

•  Multiple  scattering  with  cylindrical  geometry 

•  Three-dimensional,  non-uniform,  axisymmetric  medium 

•  Arbitrary  particle  mixtures  and  size  distributions 

•  Arbitrary  angular  scattering  phase  function 

•  Line-by-line  and  band-model  treatments 

•  Alternate  levels  of  engineering  approximation 

The  approach  to  meet  these  objectives  is  described  in  this  section. 
Our  efforts  included  (i)  the  selection  of  physical  models  and 
solution  methods  for  the  gas/particle  radiation  transfer  problems, 
(ii)  the  justification  of  selected  models  and  methods  by  means  of 
diagnostic  calculations,  and  (iii)  the  final  development  of  the 
selected  models  and  methods  for  integration  into  SIRRM. 

In  Section  2.2,  we  provide  the  requisite  technical  back¬ 
ground,  beginning  with  the  equation  of  radiative  transfer  and  its 
implicit  and  explicit  solution  schemes. 
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In  Section  2.3,  we  review  and  assess  the  combined  inter¬ 
disciplinary  experience  of  contributors  in  several  fields  (e.g., 
planetary  and  stellar  atmospheric  physics,  radiant  heat  transfer, 
plume  physics)  where  multiple  scattering  techniques  have  reached 
an  advanced  stage  of  development.  Based  on  this  assessment  we 
se1ect  the  most  appropriate  approach  to  the  plume  problem,  and 
Justify  this  selection  on  the  basis  of  key  physical  requirements 
and  model  capabilities.  The  discussion  presented  in  this  section 
formed  one  of  the  essential  elements  of  the  Program  Plan  (Data 
Item  1). 

The  other  essential  element  of  the  Program  Plan  was  a  set 
of  recommendations  for  the  final  development  of  the  scattering 
model  to  be  incorporated  into  SIRRM.  Since  January  25,  1979,  the 
recommendations  were  followed  and  the  results  of  developing  en¬ 
gineering  approximations  and  the  N-flux  model  are  reported  in 
Sections  2.4  and  25,  respectively,  together  with  numerical  results 
of  diagnostic  and  parametric  calculations. 

The  comparisons  of  the  computational  results  between  the 
different  models  developed  in  the  present  program  and  with  in¬ 
dependent  model  calculations  are  made  in  Section  2.6.  Conclusions 
about  the  models  and  the  results  are  summarized  in  Section  2.7. 

2 . 2  Scattering  Theory 

The  transfer  of  radiant  energy  through  matter  that  emits, 
absorbs,  and  scatters  radiation  was  first  described  by  astro¬ 
physicists  at  the  beginning  of  this  century.  The  governing 
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equation  was  derived  by  specializing  the  Boltzmann  equation  of 
kinetic  theory  to  the  case  of  photons.  To  this  day,  the  original 
formulation  remains  essentially  unchanged.  However,  the  range 
of  problems  amenable  to  exact  solution  is  quite  small,  and  for 
most  cases  of  interest  an  approximate  treatment  is  required. 

The  formal  solution  to  the  equation  of  transfer  provides  the 
essential  framework  for  development  of  a  coupled  gas/particle 
radiation  model.  The  task  is  to  obtain  a  practical  mathematical 
solution  based  upon  reasonable  physical  approximations  to  the 
spectral  and  spatial  characteristics  of  the  medium.  In  this 
section  we  review  the  available  computational  methods  in  order 
to  identify  promising  approaches  for  the  problem  at  hand. 

A  comprehensive  discussion  and  development  of  the  radiant 
transfer  equation  was  given  by  Viskanta  (1964) .  This  equation 
may  be  written  in  che  integro-dif ferential  form 


where 


dNx(!) 

3s 
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spectral  radiance  at  point  s  in  direction  s 
phase  function  for  scattering  from  s'  to  s 
blackbody  spectral  radiance  at  point  s 
cross  sections  for  absorption  and  scattering 
number  density 
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The  terms  on  the  right-hand  aide  represent,  respectively,  the 
extinction  (absorption  plus  scattering-out)  and  production  (emission 
plus  scattering-in)  of  radiation  in  the  direction  s.  The  phase 
function  is  normalized  such  that 
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This  formal  solution  is  not  explicit  in  that  the  local  radiance 
is  seen  to  depend  on  the  radiance  distribution  throughout  the 
medium  which  is  not  known  at  the  outset.  Explicit  solutions 
can  be  obtained  only  for  special  cases  as  noted  below. 

In  the  absence  of  scattering  (a<s)  *  0),  the  formal 
solution  to  the  equation  of  transfer  assumes  the  explicit  form 

.  s  dT  (s '  s) 

N,(s)  =  N,  (s  )T(s  ,  s)  +  /  N^s1)  —  i - —  ds' 

A  A  o  o  A  ds ' 


Tx(s',s)  »  exp |tx(s')-tx(s) 

and  where  the  optical  depth  for  extinction  and  the  albedo  for 
single  scattering  are  defined  as 

tx(s)  *  /3n(ax(a)+cx(s))ds’ 
a  (s) 
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vhere  the  transmittance  between  points  a'  and  s  la  given  by 


The  evaluation  of  this  equation  on  a  monochromatic  or  spectrally- 
averaged  basis  is  the  end  object  of  numerous  existing  line-by- 
line  or  band  model  radiation  codes. 

The  finite  difference  form  of  the  formal  solution  to  the 
equation  of  transfer  may  be  written  directly  from  the  preceding 
integral  equations.  Thus,  the  radiance  at  position  k  in 
direction  i  is  given  by 


N 


ik  "  Ni,k-lTk-l,k  +  %(Sik  +  Si,k-1)  (1_Tk-l,k) 


(2-7) 


(i”l , 2 , . . . , Iik-1,2, . . .K) 


where  the  source  function  is  given  by 

I 

sik  *  O'-iX  +  “k  £  W 


J-l 


(i"l ,2, . . . , I ;k«l , 2 , . . . K) 


and  where  the  transmittance  and  phase  integral  are  given  by 


rk-i.k  ■  exp 


<ak-i  +  °k-i)nk-isk-i  -(oka)  +  °ks 


>)nk6k] 


(k-1,2 . K) 


>  -  i  f 

jf  4tt  P(sj,si)dfi(sj) 


(i,j-l,2 . 1) 
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In  these  equations,  the  subscripts  i.j.k  denote  the  direction 
s, s'  and  position  s,  respectively,  and  the  subscript  X  has 
been  deleted  for  clarity.  With  appropriate  boundary  conditions 


on  the  radiance  N^^(k"0 , K; i-1 , 2 . . . , I) ,  which  include  the  effects 
of  any  external  sources,  these  equations  constitute  a  closed 


Solution  techniques  may  be  grouped  into  two  fundamental 
classes:  implicit  or  explicit.  The  implicit  methods  solve  the 
coupled  equations  for  and  Sik  in  an  iterative  manner 
starting  from  an  assumed  initial  condition.  A  convenient  proce¬ 
dure  is  to  define  an  initial  source  function  equal  to  the 

known  thermal  emission  component  and  corrected  source 

(2) 

function  directly.  This  procedure  is  then  repeated  to 

convergence  according  to  the  iteration  formulae 


(n)  (n)  (n)  (n) 

Nlk  -  Nl.k-1  Tk-l,k  +  *<Sik  +  St,k-1>  (1*Tk-1,k) 


(n+1) 


P.  .N.. 
Ji  Jk 


i-1.2,. 
k-1,2,. 
'n-1,2, . 


.  ,1 

• >  K  (2-8) 

.  ,N 


where  the  superscript  (n)  denotes  the  number  of  iterations  or, 
equivalently,  the  order  of  scattering. 
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The  explicit  methods  solve  a  single  matrix  equation 
for  obtained  by  direct  substitution  of  the  equation  for 
S^.  After  rearrangement  of  terms,  the  equation  may  be  written 
in  the  form 


j^|Njk<25jiT  k-l,k*“kPji* 

+  a-v-iK-i 


Nj  ,k~l(26ji  +  “k-l^i5  " 


(1=1,2,.  .  .  ,I;k-l,2 . K) 


(2-9) 


where  <5^  is  the  Kronecker  delta. 

Implicit  and  explicit  methods  offer  different  relative 
advantages  depending  upon  the  optical  and  geometric  properties 
of  the  medium.  These  properties  control  the  number  of 
iteration  steps  (implicit  methods)  or  the  size  of  the  inversion 
matrix  (explicit  methods)  required  to  obtain  a  given  level  of 
accuracy. 

2 . 3  Review,  Assessment  and  Selection  of  Multiple 
Scattering  Formulation 

2.3.1  Review  and  assessment.  In  the  presence  of 
scattering,  explicit  solutions  can  be  obtained  only  for 
suitable  approximations  to  the  radiation  flux  or  to  the 


phase  function.  Approximations  to  the  radiation  flux  in¬ 
clude  the  optically  thin  limit  (in  which  the  medium  is  not 
dense  enough  to  absorb  or  scatter  radiant  energy  from  other 
elements  in  the  volume)  and  the  optically  thick  limit  (in  which 
the  medium  is  so  dense  that  radiant  energy  from  other  elements 
does  not  penetrate).  These  limits  are  of  interest  only  as  test 
cases  for  a  more  general  model  of  arbitrary  optical  opacity. 
Apprroximations  to  the  phase  function  include  the  discrete- 
ordinate  method,  the  moment  method,  and  the  method  of  expansion 
in  orthogonal  functions  (Krook,  1955)  from  which  the  well-known 
two-stream  approximations  of  Schuster-Schwarzschild  and  Milne- 
Eddington  (Viskanta,  1964)  may  be  derived  as  special  cases. 

In  addition  to  these  explicit  methods,  various  implicit  solution 
techniques  have  been  developed.  These  techniques  typically 
construct  a  convergent  solution  by  incremental  steps  in  space  or 
time. 

A  summary  of  available  multiple  scattering  solution  techniques 
is  given  in  Table  2-1.  These  methods  can  be  roughly  grouped 
into  four  general  classes:  superposition,  transformation, 
convergence,  and  matrix.  Each  of  these  is  discussed  and  developed 
in  detail  in  the  survey  papers  by  Hansen  and  Travis  (1974)  and 
Irvine  (1975).  The  identifications  are  necessarily  somewhat 
arbitrary  in  that  some  interrelated  techniques  are  included 
and  some  distinguishable  subsets  are  excluded.  Tables  2-2  and 
2-3  list  selected  key  investigations  according  to  solution 
technique  and  chronology,  respectively.  The  former  was  compiled 
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Addition 
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Oyntltaaia 


X  and  Y 
functions 


Spherical 
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TABLE  2-1.  METHODS  OF  SOLUTION  FOR  MULTIPLE  SCATTERING 


1.  Superposition  Methods 


Sups  rpoeit ion  of 
layers  of  known 
transmission  end 
reflection 

Superposition  of 
optically  thin 

layers 

Superposition  of 
volume  elements 


Wide  application 
for  planetary 
atmospheres 


lior.  -  homogeneous 
atmospheres 


Generality 


Mathematical 
transformation 
of  variant 
imbedding 

Phase  integral 
replaced  by  »u u. 
of  sphurical 
harmonics 


Transformation  Method* 

.'Single  angle 
dependence 


General 

applicability 


Mcn-aultting 
medium i  plsisa- 
parellel  slab 


Long  computations 


Inefficient  for 
syammtric  media 


Rayleigh 

scattering 


Large  matrix 


Eigenfunction 


Phase  lntagrai 
replaced  by 
sum  of  orthogonal 
eigenfunctions 


Ha  themat  i  cal  ly 
rigorous 


Numerically 

prohibitive 


3.  Convergence  Methods 


Iteration 


Successive  down¬ 
ward  and  upward 
integrations  to 
convergence 


Simple  concept 


Low  edacity 


Successive 

orders 


Superposition  of 
solutions  for 
different  orders 
of  scattering 


Simple  concept 


Ltw  opacity; 
non- emitting 
medium 


Monte  Carlo 


Individual  photon 
histories  traced 
end  superposed 


Sisf>le  concept; 
general 
•pplicabil ity 


long  computation* 


<•  Matrix  Methods 


Discrete 

ordinate 


Phase  Integral 
replaced  by 
finite  sum 


General 
eppllceb.il  ity 


Large  matrix 


finite 

difference 


operates  an 

differential 

equation 


Ooap  rational 
experience 


Large  matrix 
unless  sysmmtric 


finite 

eleemnt 


Gperetee  on 
integral  equetlo 


Reduces  to 
two-stream  Umlt 


Large  matrix 
uftieaa  nyemmtric 
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TABLE  2*2,  REFERENCES  FOR  MULTIPLE  SCATTERING  METHODS 


METHOD 


REFERENCES 


Layer  Addition 


van  de  Hulst,  et  al.  (1963,  1968)  ,  Twomey, 
et  al.  (1966),  Hausen  (1969,  1971),  Hovenier 
(1971),  L&cis  and  Hansen  (1974),  Grant  and  Hunt 
(1969,  1971),  Plans,  et  al.  1.1973), 
Preisendorfer  (1965) . 


Successive  Orders 


van  de  Hulst  and  Irvine  (1962) ,  Irvine  (1964, 
1965) ,  Uesugi  and  Irvine  (1970) ,  Poon  and  Ueno 
(1974),  Dave  (1564) ,  Hovenier  (1971),  Chou 
(1978) . 


Invariant  In-bedding 


Iteration 


Uesugi  and  Irvine  (1970),  bellnun,  et  al.  (1960, 
1963,  1967,  1969),  Wing  (1562). 

Herman,  et  al.  (1963,  1965,  1970),  Dave,  et  al. 
(1970),  Crcsbie  and  Linsenbardt  (1975). 


X  and  Y  Functions 


van  de  Hulst  (1970) ,  Sobolev  (1972,  1974) , 
Carlstedt  and  Muilikin  (1966) ,  Pahor  and  Kuscer 
(1966),  Busfcridge  (1960,  1967). 


Discrete  Ordinate 


Chandrasekhar  (1950),  Jefferies  (1955),  Chu  and 
Churchill  (1955) ,  Lenoble  (1956) ,  Chin  and 
Churchill  (1965),  Samuelson  (1969),  Liou  (1973, 
1974),  Weinman  and  Guetter  (1972),  Whitney  (1972, 
1974),  Kofink.  (1967),  Love  and  Grosh  (1965). 


Spherical  Harmonics 


Canosa  and  Penafiel  (1973) ,  Davison  (1958) ,  Case 
and  Zweifel  (1S67) ,  Deuze,  et  al.  (1973),  Dave, 
et  al.  (1974). 


Eigenfunctions 


Case  and  Zweifel  (1967) ,  Kuscer  and  McCormick 
(1973,  1974),  Kaper,  et  al.  (1970),  Mika  (1961). 


Monte  Carlo 


Plass  and  Kattawar  (1968,  1971),  Danielson, 
et  al  (1968,  1969) ,  Sanford  and  Pauls  (1973) , 
Collins,  et  al.  (1965,  1972),  Van  Blerkom,  et  al. 
(1.971,  1974),  McKee  and  Cox  (1974),  Kattawar, 
et  al.  (1971,  1973),  Marchuk  and  Mikhailov  (1967) 


Dj  f fusion 


Heasley  (1977) ,  Barkstroo  and  Arduini  (1978) , 
Laung  (1975,  1976) 


Finite  Difference 


Barkstrom  (1976) 
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TABLE  2-3 .  CHRONOLOGY  OF  MULTIPLE  SCATTERING  METHODS 


YEA*  PEfEHEKCE 


PHYSICAL  SITUATION 


KFTHOD  OF  ANALYSIS 


1955  Jefferies 


Point  and  Infinite  lint  sources 
in  a  sanl-lnflnlte  nedlun 


Diffusion  approximation 


1355  Jifftriis 


1955  I  Chu  me 
I  Churchill 


A  atrip  of  collimated  radiation 
lreldent  on  a  senl-lnflnlta 
■ad lum 

Foil;**,  aourc*  outiide  an  sn'so- 
tropically  scitterl nr  finite 
layer 


Discrete  ordinates  and 
Fourier  transform 


Six-flux  approximation 


1956  I  Qlovanelll 


An  infinite  line  source  on  the 
surface  or  within  a  coml-lnrin- 
lte  medium  and  a  thick  finite 
layer 


Approach  of  Jefferies 


1358  I  Chandrasekhar 


An  Infinitely  narrow  eolllnated 
beam  Incident  on  a  senl-lr finite 
medium 


Invariance  principles 


l«5l  Bellman,  Kalaba, 
and  Uen© 


An  Infinitely  narrow  collimated 
beam  incident  on  a  finite  lever 
with  oronerttes  varying  with 
depth 


Invariant  Imbedding 


1963  I  Bellman,  Kalabs, 


196*  I  Smith 


Collirafced  radlntlor  incident 
on  «n  inhomogeneous  unlrotmrl- 
cally  scattering  finite 
1 

A  unlferm  strip  of  collimated 
radiation  incident  on  a  seml- 
inflntte  medium 


Invariant  Imbedding 


Assured  cosine-varying  In¬ 
cident  radiation  to  reduce 
the  two-dimensional  In¬ 
tegral  equation  for  the 
source  function  to  a  one- 
dlmenalonal  form 


196*  j  Nalkevlch,  Fonln, 
J  and  Hoxenherg 


Collimated  radiation  Incident  on  Horliontal  variation  of 
a  finite  layer  with  nonunlforat  downward  radiation 
reflection  from  lowvr  surface  neglected 


1965  I  CMn  and 

ft  Churchill 


An  arbitrary,  cylindrical  source  Six-flux  approximation 
in  ar  ani sotroplcal ly  scattering 
finite  slab 


1966  I  Dolln 


Harrow  collimated  beam  Incident 
on  a  seml-lnfinlte  anlaotropl- 
cally  scattering  medium 


Small  angle  transport  equa¬ 
tion  used  for  unscsttered 
radiation,  diffusion  equa¬ 
tion  used  for  scattered 
radiation 


1967  R  Vtlnman  and 
I  Swart trauber 


Collimated  radiation  Incident 
on  a  nonabaorMng  finite 
layer  with  horlxontally  non- 
uniform  scattering  coefficient 


Qlovanelll (a  approximation 
used  to  solve  transport 
equation 


1967  I  Smith  and  Hunt 


1967  I  Hunt 


A  uniform  strip  of  collimated 
radiation  Inclder.t  on  a  senl- 
flnlte  layer 

Collimated  radiation  incident 
normally  on  an  anlaotroplcally 
scattering  finite  redlum 


Same  at  Smith  196* 


Expression  or  scattering 
function  as  series  of 
Legendre  polynomials 


1967  I  ftobeo 


An  Isothermal  ueml-xnfinltt  slab  Diffusion  equation 

solved  by  iterative 
technique 


19(7  I  Drobyshevlch 


Collimated  raulatlon  Incident  on  Horizontal  variation  of 
an  anlsotroplcally  scattering  downward  radiation 

finite  layer  with  nonuniform  neglected 

reflection  from  lover  surface 


1960  I  Pomanova 


An  Infinitely  narrow  collimated 
beam  normally  Incident  on  s 
•ami-infinite  anlaotrooi rally 
scattering  medium 


Three-dimensional  moments 
used  to  solve  the  trans¬ 
port  equation  for  large 
depths 


I960  j  Pomsnova 


An  irifln4telv  narrow  colligated 
bear  normally  Incident  on  an  in- 
Isotropically  scattering  finite 
lever 


Frail  angle  aoproxlmat Ion 
to  the  transnort  equation 
solved  using  spatial 
rorents 
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TABLE  2-3.  (CONTINUED) 


mi 

rtmutt 

physical  ijtvatio* 

rmtoD  or  a-altiu 

i»« 

rank 

Collimated  eoi  nt  source  norrally 
incident  on  an  ar.l  aotropleali* 
scattering  finite  layer 

Font*  Carle  teetoilqiN 

i«i 

Hunt 

Collimated  Pease]  function 
varying  radiation  incident 
on  a  finite  layer 

Internal  eauatlon  for  the 
source  function  formulated 
using  Green’s  function 

!•<» 

FtoeFM"  and 

Lott 

An  Isothermal  finite  cylindrical 
medium;  both  laotroflc  ard  ani¬ 
sotropic  scatterirr 

vonte  Carlo  technioue 

'.»<* 

Hunt 

Colltrafced  radiation  rorriHv 
Incident  on  a  finite  layer 

Integral  equation  for 
source  function  trans¬ 
formed  Into  singular 
eouatlon  for  the  re¬ 
solvent  kernel  using 
properties  of  a  Green’s 
funct ion 

1*«? 

Brave-ZMvotov.ikiy , 
Dclln,  Luchlnln, 
and  £*r»rye» 

Narrow  collimated  bear  norrally 
Incident  on  a  aenl-lnflnlte  an- 
isetroplcallt  reattermr  rellur 

Uses  solution  to  transport 
equation  found  In  Do! In  19 66 

1S«* 

GolubltsMy  and 
Tantosnev 

An  Infinitely  narrow  collirated 
lean  at  the  center  of  an  anlso- 
trorlcally  ccattrrlnr  srherical 
red  1  urn 

yionte  Carlo  technique 

1  *J70 

Love  and  Turner 

An  l»oth*rral  serl-inf lrlte  slab 

Ponte  Carlo  technique 

IfTO 

Scofield  end 

Love 

A  atrip  of  dlffus**  radiation  In¬ 
cident  on  an  **rlttlnr  aerl-ln- 

flnlte  medium,  both  Isotropto 
and  anisotropic  scstterlnr 

rente  Carlo  technioue 

1*71 

Hjlllcltl 

An  infinitely  narrow  collimated 
beam  incident  on  a  finite  layer 
and  ser !-inf lnlt*  medium 

Fourier  transform  and  In¬ 
variance  technloues 

1*71 

*on»rovm 

An  Infinitely  narrow  coll  trotted 
bean  normally  Incident  on  a  non- 
absorbing  rinlte  layer 

Three-dimensional  Invariance 
principle®  used  to  obtain 
equations  for  spatial  mo¬ 
ment  c  of  intensity 

1*71 

Homonova ' 

An  infinitely  narrow  colllrated 
bean  normally  Incident  on  an  an- 
* sotroptcelly  acatt-'ring  finite 
layer 

Spatial  moments  used  to 
solv'-  ana  11  arifle  apcroxl- 
natlon  to  the  transport 
enuat Ion 

1971 

Koeht tkpv 

An  Infinitely  narrow  collimated 
beo*>  normally  incident  on  *r>  an- 
laotroplcally  scattering  finite 
layer 

method  of  spherical  harmonica 

1*77 

Bushnakova,  Z«*«, 

and  Kataev 

Diffuse  and  unidirectional 
source*  inside  an  ar.lsotrcplcal- 
ly  scattering  finite  layer  and 
seml-lnf lnlte  medium 

Hankel  transformation  used 
to  solve  diffusion  equation 

1*7  3 

Fomnnova 

An  Infinitely  narrow  collimated 
beam  normally  Incident  on  a  fi¬ 
nite  layer  with  weak  absorption 

Intensity  represented  by  trans¬ 
verse  momenta  and  invariance 
principles  used 

1*7J 

Funk 

A  narrow  collimated  bear*  and 
noncolllmated  source  in  an  In¬ 
finite  enisotroplcally  scatter¬ 
ing  medium 

Monte  Carlo  technique 

1*7J 

Cordon  and 

Bnlttel 

An  infinitely  narrow  collimated 
bear  nomallv  Incident  on  an  an- 
lsotroolcally  scattering  finite 
layer 

Iterative  ray  tracing,  no  re¬ 
turn  backrcetterlng  and  small 
angle  icpro* imatlona  used 

1*7* 

Crusble  and 

Koewing 

A  step  and  finite  atrip  of  col¬ 
limated  and  diffuse  radiation 
Incident  on  a  finite  layer 

Assumption  of  cosine  vary¬ 
ing  Incident  radiation  re¬ 
duces  problem  to  one-dl- 
stnsiona]  problem; 

AmhirrurMane  method 

1*7* 

Beckett,  Foster, 

Hutson  and  Host 

A  collimated  cylindrical  beam 
normally  incident  on  a  finite 
layer 

Hunt’s  Integral  eauatlors 
solved  by  Iteration 
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from  previous  review  papers  together  with  our  own  review; 

the  latter  was  taken  directly  from  Crosbie  and  Llnsenbardt 

(1976).  Neither  list  is  comprehensive,  but  both  serve  to 

illustrate  the  depth  and  scope  of  previous  work.  Most  of 

these  studies  applied  to  non-emitting  plane-paralJ.el  atmospheres. 

Only  a  limited  number  of  the  studies  included  the  effects  of 

* 

thermal  emission,  inhomogeneous  media,  non-isotropic  scattering, 
and  non-planar  geometry  which  are  important  in  the  present 
application.  None  of  the  studies  treated  all  these  effects 
simultaneously.  Despite  these  deficiencies,  this  previous  work 
provides  a  firm  base  for  the  selection  and  development  of  a 
multiple  scattering  model  for  two-phase  plumes. 

Existing  rocket  plume  radiation  codes  which  treat  particle 
radiation  were  also  reviewed  and  are  presented  in  Table  2.4, 

The  early  codes  of  Fontenot  (1965),  Gulrajani  (1964),  Hunt  (1966) 
and  the  Aeronutronic  group  (Carson,  1965)  were  all  developed  to 
predict  base  heating.  The  more  recent  plume  codes  of  Wilson 
(1973) ,  Rieger  (1974) ,  Vanderbilt  and  Slack  (1976) ,  and  the  ARAP 
group  (Fishburne,  1977)  were  developed  to  predict  plume  IR  signatures 
These  codes  approximate  particle  radiation  effects  by  either  of 
two  methods:  Pseudo-gas  or  one -dimensional  beam.  The  pseudo-gas 
approximation  includes  particle  emission/absorption  effects  but  negle 
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TABLE  2-4.  EXISTING  ROCKET  PLUME  RADIATION 
CODES  INCLUDING  PARTICLES 


PLUME  CODES 

METHOD 

COMMENTS 

Fontenot 

(1965) 

Inverse  wavelength 
method 

Simplistic  approach;  single  particle 
size  neglects  scattering.  Estimates 
base-heating  for  constant  radial 
properties . 

Gulrajani  (1964) 
Hunt  (1966) 

Neutron- scattering 
analogy  method 

Based  on  Morizumi  and  Carpenters 
(64)  analysis.  Five  particle  sizes. 
Includes  1-D  multiple  scattering. 
Neglects  anisotropic  scattering 
effects.  Assumes  Qt-*-2.  Base-heating 
estimates . 

Aeronutronic 
(Carlson,  1965) 

One-dimensional  beam 
approximation 

Plume  codes  to  analyze  Saturn  ullage 
retro  and  strap-on  exhaust  plumes. 

1-D  anisotropic  scattering  approxi¬ 
mation. 

Wilson  and  Hahn 
(1973) 

One-dimensional  beam 
approximation 

Uncoupled  gas  and  particles.  Adds 
separately  computed  particle  radia¬ 
tion  and  gas  radiation.  Narrow 
inband  resolution.  Treats  5  particle 
sizes . 

Rieger  (1974) 

Pseudo-gas  approxi¬ 
mation 

Computes  coupled  particle  and  gas 
emission  and  absorption. 

Vanderbilt  & 

Slack  (1976) 

One-dimensional  beam 
approximation 

Couples  multiple  scattering  and 
gaseous  emission/absorption  along 
beams.  Spectral  resolution  5  cm-1. 
Treats  5  particle  sizes.  Treats 
gases  as  non-scattering  particles 
yet  retains  correct  solution  of 
radiative  transfer  equation  for  gases 
alone  and  particles  alone.  Spectral 
and  inband  output  coupled  to  LOWTRAN 
codes.  nave  code  for  Q  ,  Q 

S  cl 

ROCRD  Code 
(ARAP) 

Pseudo  -gas 
approximate  n 

Includes  particle  emission,  but  neg¬ 
lects  scattering.  Assumes  Beer's 
law. 
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scattering.  The  one -dimensional -beam  approximation  includes 
particle  scattering  effects  along  the  beam  but  neglects  scattering 
into  the  beam  from  other  directions.  Of  the  two  methods, 
the  one -dimensional  beam  provides  a  better  physical  description 
of  particle  radiation  effects  with  no  appreciable  penalty  in 
computational  complexity.  Wilson  used  the  one- dimensional  beam 
approximation  to  compute  the  particulate  radiation  field  in  an 
uncoupled  manner  for  subsequent  addition  to  the  gaseous  radiation 
field;  the  validity  of  this  procedure  is  restricted  to  low 
optical  depchs.  Vanderbilt  and  Slack  used  the  one - dimens ional 
beam  approximation  to  compute  the  coupled  particulate  and  gaseous 
radiation  fields,  thereby  allowing  the  gas  to  absorb  the 
particulate  continuum  and  the  particles  to  scatter  and  absorb 
the  gaseous  radiation.  The  code  developed  by  Vanderbilt  and 
Slack  is  the  best  of  the  existing  codes  for  predicting  spectral 
IR  signatures  from  missile  plumes.  However,  the  lack  of  coupling 
between  adjacent  1-D  beams  means  that  it  does  not  meet  the  SIRRM 
code  requirements,  in  general. 

One  dimensional  beam  methods  have  received  extensive 
application  and  limited  verification  .in  the  area  of  radiative 
heat  transfer  but  only  limited  application  and  no  verification 
in  the  area  of  radiative  signatures.  This  distinction  is 
important  because  of  a  fundamental  difference  between  the  two 
phenomena.  On  the  one  hand,  the  radiative  heat  transfer  to 
an  adjacent  point  is  the  sum  of  contributions  from  non-parallel 
beams.  On  the  other  hand,  the  radiative  signature  to  a  remote 
point  is  the  sum  of  contributions  from  parallel  beams.  The 
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latter  effect  is  much  more  sensitive  to  errors  in  the  directional 
properties  of  the  radiation  field.  Thus,  the  SIRRM  code  must 
include  a  three-dimensional  scattering  model  in  order  to  achieve 
the  required  directional  accuracies. 

2.3.2  Selection  and  justification.  As  noted  previously, 
the  SIRRM  scattering  model  for  two-phase  plumes  must  include  the 
following  key  physical  effects:  thermal  emission,  inhomogeneous 
media,  non-isotropic  scattering,  and  non-planar  geometry.  In 
addition,  the  SIRRM  scattering  model  must  possess  the  following 
computational  attributes:  accuracy,  efficiency,  and  generality. 
These  general  attributes  have  rather  specific  meanings  for  the 
present  application  as  noted  below: 

Accuracy  The  chosen  technique  should  incorporate  sufficient 
spatial  resolution  to  yield  an  accurate  mathematical  solution. 
The  technique  should  demonstrate  adequate  agreement  with 
the  predictions  of  alternative  methods  for  suitable  test 
cases . 

Generality  The  chosen  technique  should  incorporate  suf¬ 
ficient  spatial  dimensions  to  represent  the  real  physical 
problem.  The  technique  should  treat  all  ranges  and  varia¬ 
tions  of  optical  depth,  scattering  albedo,  and  phase  func¬ 
tion  such  that  no  limit  is  imposed  on  the  nature  of  gas 
or  particle  properties,  plume  properties,  or  observation 
geometry. 


Efficiency  The  chosen  technique  should  yield  computed 
results  in  machine  times  which  are  comparable  to  current 
gaseous  emission  codes  for  an  equivalent  number  of  aspect 
angles.  Approximate  techniques  which  incorporate  physically 
acceptable  assumptions  should  be  evaluated  in  order  to 
determine  the  most  efficient  technique  for  different  classes 
of  problems. 

The  available  methods  were  screened  according  to  these  criteria. 

The  dual  requirement  for  engineering  efficiency  and  scien¬ 
tific  accuracy,  together  with  the  need  for  computational  gen¬ 
erality,  and  flexibiltiy,  led  to  the  selection  of  straight¬ 
forward,  finite  difference  procedures  of  variable  spatial  and 
angular  resolution  applied  to  locally  cylindrical  plume  geometries. 
Two  basic  approaches  (herein  identified  as  "6-flux"  and  "N-flux") 
were  distinguished  and  various  alternatives  were  developed  within 
each  classification.  The  semantic  distinction  refers  to  the 
associated  subdivision  of  angular  space  (6  orthogonal  coordinates 
versus  N  non-orthogonal  coordinates  in  azimuth  and  elevation) . 
However,  the  fundamental  distinction  li^s  with  the  generality 
of  the  N-flux  model  versus  the  computational  simplicity  of  the 
6- flux  model. 

Within  the  6-  flux  model  class,  two  independent  sub-class 
models  were  developed.  The  first  model  employs  suitable  approxi¬ 
mations  for  the  axial  and  tangential  radiances  in  order  to 
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reduce  the  solution  to  a  one-dimensional  radial  integration 
similar  to  existing  two-stream  models  (Vanderbilt-Slack,  1976). 

The  second  model  employs  a  transformation  of  the  multiple¬ 
scattering  source  function  in  order  to  reduce  the  solution  to 
a  coupled  set  of  radial,  tangential,  and  axial  integrations,  each 
of  which  is  functionally  similar  to  existing  two-stream  models. 

The  two  six-flux  models  are  described  and  compared  in  Section 
2.4. 

The  N-flux  model  class  includes  a  generalized  set;  of 
optional  finite-element  grid  geometries  including  independent 
specification  of  the  angular  resolution  in  the  azimuth  and 
elevation  planes.  The  model  solves  the  coupled  three-dimensional 
equations  of  transport  using  a  generalized  matrix  network 
to  describe  the  radiance  distribution  and  source  function  at:  each 
grid  point.  The  model  can  be  exercised  as  either  an  engineering 
or  scientific  tool  by  simple,  variations  in  the  spatial  and  angular 
grid  network  to  achieve  low  or  high  resolution.  Analytical 
details  are  presented  in  Section  2.b. 

Each  of  the  models  was  exercised  parametrically  over  the 
full  physical  range  to  illustrate  the  effect  of  optical  depth, 
scattering  albedo,  and  phase  function  on  the  local  and  integrated 
target  signature  characteristics.  These  calculations  provided 
independent  tests  of  the  physical  credibility  and  computational 
efficiency  of  the  separate  models  as  well  as  comparative  tests 
of  the  accuracy  and.  se  Lf-consistency  of  the  composite  model  set. 
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Additional  verification  was  accomplished  by  performing  inde¬ 
pendent  calculations  with  the  NASA  plane-parallel  slab  code 
(Barkstrom,  1976)  modified  to  include  self -emission  and  with 
the  Grumman  Monte  Carlo  code.  Results  of  these  calculations 
are  summarized  in  Section  2.6. 


Six-Flux  Ent 


The  basis  for  the  engineering  approaches  developed  under 
Task  IA  is  the  approximation  of  the  multiple-scattering  source 
function  by  a  six-element  orthogonal  angular  quadrature.  In 
this  approximation,  the  general  differential  equation  of  transport 


-  -N(T,ft)+[l-w(T)|  N°(x)+  yP(T,fl,n')N(T  ffl')dr<'  (2-10) 


is  written  for  each  orthogonal  direction  (i  =  1,2,3)  in  the  specific 
form 
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where  the  terms  on  the  right-hand  side  represent  the  respective 
contributions  of  extinction,  emission,  and  multiple  scattering 
to  the  radiance  gradient  and  where  the  terms  f,b,s  represent 
the  forward,  backward,  and  sideward  components  of  the  phase  func¬ 
tion  (f+b+4s»l)  This  quadrature  provides  the  minimum  angular 
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resolution  necessary  to  account  for  three-dimensional  multiple¬ 
scattering  effects.  For  comparison  to  previous  approaches,  the 
standard  two-element,  one-dimensional  approximation  is  obtained 
as  a  special  case  (s=0) . 

The  orthogonal  coordinates  are  chosen  to  coincide  with  the 
natural  cylindrical  coordinates  (r,8,z)  of  the  axisymmetric 
plume.  Figure  2-1  illustrates  the  grid  network  formed  by  the 
intersections  of  these  coordinates  with  the  boundaries  of  the 
finite  elements  (radial  and  axial) .  The  six- flux  equations  are 
solved  as  a  coupled  set  along  these  mutually-perpendicular 
integration  paths.  Radiances  for  oblique  lines-of-sight  are 
computed  in  a  subsequent  step  based  upon  the  known  six-flux 
solution. 

Within  the  six-flux  model  class,  two  different  levels  of 
engineering  approximation  were  employed.  The  first  method 
yields  an  approximate  explicit  solution  based  upon  elimination 
of  the  tangential  (9)  equations  of  transport  in  favor  of  simple 
closure  relations.  The  second  method  yields  an  exact  implicit 
solution  based  upon  iteration  or  inversion  of  the  fully-coupled 
multi-dimensional  equation  set.  Both  methods  were  originally 
developed  for  infinite  cylindrical  geometries  (negligible  axial 
gradients) .  The  second  method  was  subsequently  extended  to 
finite  cylindrical  geometries  in  order  to  assess  the  influence 
of  plume  axial  gradients  on  the  radiance  at  side  aspect.  The 
following  paragraphs  provide  analytical  details. 
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2.4.1  Method  I.  This  method  is  based  upon  a  transformation 
of  the  three-dimensional  equations  of  radiative  transport  to  equi  - 
valent  one -dimensional  form  for  infinite  planar  or  cylindrical  media. 
The  resulting  equation  is  integrated  along  a  piecewise-uniform 
path  perpendicular  to  the  surface  to  yield  an  explicit  solution 
for  the  normal  and  tangential  radiances  throughout  the  medium. 
Radiances  at  intermediate  angles  are  determined  by  analytical 
interpolation.  No  restriction  is  placed  upon  the  range  of 
physical  properties  considered  (e.g.  optical  depth,  scattering 
albedo,  phase  function)  or  upon  the  variation  of  those  properties 
within  the  medium. 

The  ability  to  treat  planar  and  cylindrical  geometries  in  a 
self-consistent  manner  provides  an  important  connection  between 
the  available  literature  (mostly  planar)  and  the  plume  problem 
(nearly  cylindrical).  As  a  result  of  this  connection,  the  utility 
of  existing  plane  parallel  codes  (e.g.,  Barkstrom,  1976)  for 
defining  baseline  parametric  trends  is  greatly  enhanced.  In 
addition,  solutions  for  planar  media  which  are  infinite  in  one 
or  two  directions  (i.e.,  ribbon  or  slab,  respectively)  provide 
convenient  bounds  on  the  emitted  radiation  from  equivalent 
cylindrical  media. 

In  the  rectangular  r,0,z  coordinate  system  of  Figure  2-1,  the 
equation  of  transfer  reduces  to  the  six  simultaneous  equations 
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These  six  equations  are  further  reduced  to  two  (±r)  by  means  of 
suitable  analytical  approximations  for  the  perpendicular  radiances 
(±8,±z).  The  z-equation  is  simplified  through  the  assumption  of 
negligible  axial  gradients 

dN,;4. 
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z 

which  reduces  the  equation  to  simple  algebraic  form.  The  0 -equation 
is  eliminated  through  the  assumption  of  approximate  proportionality 
factors 


N0±  ~  (plane) 

N0~  «  ^f,%(Nj:++Nr"’)  (cylinder) 

which  are  piecewise-cons cant  over  a  radial  integration  step.  The 
range  0<$p<l  defines  the  limits  of  thin  and  thick  optical  paths 
in  the  8-direction.  In  general,  the  factor  $>  is  approximated 
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by  the  corresponding  radiance  ratio  for  a  non-scattering  medium. 
For  a  uniform  medium,  this  ratio  is  given  simply  by  the  emissivity 
ratio 


1-exp (-Tq) 

[l-exp(~T0)j  /  [l-exp(-t) 


where  r  is  the  total  optical  depth  for  a  cylinder  of  diameter  2R 

and  t q  is  the  half  optical  depth  measured  perpendicular  to  the 

2  % 

line-of-sight .  For  a  uniform  cylinder  r0  -  T[].-(r/R)  ]  . 

Substitution  of  these  approximations  into  the  axial  equation  of 
transport  yields  the  explicit  formulae 
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With  these  approximations,  the  six- flux  radial  equation  of 
transport,  may  be  written  in  the  functional  form  of  the  two-flux 
equation  with  modi. fled  coefficients.  That  is, 

dN  .  '  +  . 

±  “  "«Nr“  +  BN.r+  +  yN°  (2-14) 

where  the  piecewise-constant  coefficients  are  given  by 
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planar 


cylindrical 


a  =•  l-io(f+sx)  a 

0  =  a»(b+sx)  £ 

Y  =  (l+x)(l-a))  Y 

X  =  2suj(l+i'p)[l-ti,'(f+b+2<I>pS)]  ^  x 


l-(jj{f+s[x+$c(l+x)  ]} 
w{b+s  [  x+$c  (1+x)  ]  } 
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Integration  over  the  uniform  region  from  station  (k-1)  to  k  yields 
explicit  expressim\s  for  the  outward  and  inward  radiance  components 
at  the  outer  and  inner  surfaces,  respectively. 


K  -  v;  *vt-i  +%K 

N  k-l  =  CkNk  4'  plcN  k-1  +  ekNk 


(2-15) 


where  the  subscript  r  has  been  deleted  for  clarity. 

These  equations  contain  the  boundary  conditions  (N  ^  and 
^k-l>’  local  source  function  (H°)  ,  and  the  local  transport 

coefficients  (reflectivity,  transmissivity,  and  emissivity)  defined 
by  the  relations 
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The  axial  radiances  at  the  outer  and  inner  surfaces  then  follow 
directly . 

The  formulation  for  a  uniform  planar  medium  is  complete  at 
this  point.  For  a  non-uniform  or  cylindrical  medium,  the  formu¬ 
lation  is  completed  by  superposition  of  elemental  solutions  for 
K  piecewise-uniform  subregions  (k=l , 2 , . . . ,K)  subject  to  the 
overall  boundary  conditions 

=0  (no  external  sources) 

Nq  «=  Nq  (axial  symmetry) 

The  radiance  emerging  from  the  K-th  (outer)  element  is  determined 
from  Eq.  2-15  by  successive  Gaussian  elimination  beginning  at  the 
center  (k  -  0) .  The  general  recursion  relation  (developed  by 
induction)  is  given  by 
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The  surface  radiance 

(K^  =  B^)  is  determined 

explicitly  from  the 

sequential  summation. 
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The  preceding  analysis  yields  the  normal  and  tangential 
radiance  components  for  planar  and  cylindrical  media.  Between 
these  limits,  the  variation  of  radiance  with  elevation  angle  (<j>) 
or  azimuth  angle  (9)  (standard  spherical  coordinates)  must  be 
determined  by  a  subsequent  calculation.  The  following  paragraphs 
describe  a  simple  approximate  procedure  based  on  analytical 
interpolation.  An  alternative  exact  procedure  based  on  the  known 
six-flux  source  function  distribution  is  described  under  Method  It 

The  approximate  interpolation  formulae  are  based  on  weighted 
sums  of  the  local- to-normal  and  local- to- tangential  radiance 
ratios  for  a  non-scatteru.ng  medium.  From  Figure  2-2 
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— xil-  N(tt/2,tt/2)  sin% 

N<t,/2,it/2>Jw=0 


N(tt/2 , 0)  cos^<f> 

GJ=0 

where  N(Tr/2,7r/2)  and  N(tt/2,0)  are  the  normal  and  tangential 
radiances  predicted  by  the  six-flux  model.  The  0  and  <p 
variations  are  conveniently  separated  by  rewriting  N(e,<f>)  as 
the  product 
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Figure  2-2.  Estimation  of  Three-Dimensional  Radiance  Distribution 
from  Computed  Radiances  at  Normal  and  Tangent  Aspect. 
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This  separation  facilitates  the  integration  over  9  to  obtain  the 
local  station  radiation. 
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For  uniform  media  the  radiances  for  w=0  are  simply  proportional 
to  the  corresponding  emissivities ;  i.e.. 


N(9,<J)  ^  e  (rsin6/sin4>) 


where  t  is  the  optical  depth  based  on  the  cylinder  diameter. 
In  this  case,  the  radiance  is  given  by  the  simplified  formula 

N(0,4>)  =  e(Tsin0/sin4))|  N(Tr/2,Tr/2)sin20/e(T) 

+  N  (  7T  /  2 , 0)  cos^  I 


2-30 


puwuiimw* 


icWijurc" 


and  the  station  radiation  is  well-approximated  by  the  empirical 
formula  (see  Figure  2-3) 

J'(<f>)  =  2.Rsin<f>£  (i /sin<j>)  [N(tt/2  ,  tt/2)  sin^<|>/e  (t)-I-N(tt/2 , 0)cos*"<}>  ]  • 
e  (  .  546  c/sin<j>)  /  e  ( .  695t/sin4)) 

The  non-dimensional  ratio  J'  (<ji) /2Rsin<|>N(TT/2 ,  <f>)  is  observed  to 
approach  ir/4  or  1  in  the  optically  thin  or  thick  limits, 
respectively.  These  interpolatioi  formulae  suffice  to  define 
the  entire  radiation  field  at  a  given  axial  station  in  terms 
of  the  computed  normal  and  tangential  radiances. 


Figure  2-3.  Estimation  of  Station  Radiation  for  Uniform 
Cylinder  from  Computer  Centerline  Radiance . 
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2.4.2  Method  II.  This  method  is  based  upon  a  transformation 
of  the  three-dimensional  equations  of  transport  tc  equivalent  iso¬ 
tropic  form  for  infinite  and  finite  cylindrical  media.  The  resulting 
equations  are  integrated  along  and  perpendicular  to  the  radial  direc¬ 
tion  to  yield  an  implicit  solution  for  the  isotropic  (scalar)  source 
function  distribution.  Numerical  results  are  obtained  by  iteration 
using  the  method  of  successive  scatterings.  Solutions  for  arbitrary 
aspect  are  then  computed  directly  in  terms  of  the  known  source  func¬ 
tion  by  means  of  a  two- stream  integration  along  the  observer  line- 
of-sight.  As  for  Method  I,  no  restriction  is  placed  upon  the  range 
or  variation  of  physical  properties  within  the  medium. 

The  general  six-flux  equation  of  transport  (see  Eq .  2-11) 


-N.'  -f  (1-u))N°+gj  f  N  .  4-  bN/  +s  V(N.,  +  N  ",  )  (i-1, 2, 3) 

i'/i 


is  transformed  by  simple  rearrangement  of  terms  to  equivalent, 
isotropic  form 

dN  *  + 

±  =  _aN.  +  3Nt+ +  y(N°+N’)  (i-1, 2 ,3)  (2-1 6) 

where  the  piecewise-constant  coefficients  are  given  by 

a  =  l-u(f-s) 
g  =  w(b-s) 
y  -  (1-w) 

and  where  the  source  function  for  multiple  scattering  i3  given  by 
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The  key  feature  of  the  transformed  equation  set  is  that  the  source 
function  now  contains  only  the  side-scattering  component  of  the 
phase  function  such  that  N’  is  a  scalar  quantity.  This  important 
simplification  is  obtained  for  any  arbitrary  phase  function  which 
can  be  approximated  by  a  six-angle  orthogonal  quadrature. 

For  a  known  source  function  N’  (as  from  a  previous  iteration) 
the  transformed  equation  of  transport  is  functionally  equivalent  to 
Method  I  and  can  be  integrated  as  before  to  yield 


Nik  -  pikNik+  tikNik-i  +  eik<»VN,k>  ) 

(i=l , 2 , 3)>  (2-17) 

Nik-1  =  tikNik  +  pikNik-l  4  Eik(N°k+NV  ) 


where  p.^,  t^^,  are  the  reflectivity,  transmissivity,  and 

-h 

emissivity  of  the  element  k  in  the  direction  i  and  NT^, 
arc  the  surface  radiances.  The  coefficients  p,  t,  e  are 
evaluated  as  before  (see  Method  I)  in  terms  of  the  coefficients 
a,  3,  Y.  where  the  latter  have  new  definitions  as  given  above. 

The  requisite  geometrical  relations  are  given  in  Section  2.5. 

For  an  array  of  K  elements  (k=l ,  2,...K),  the  boundary  conditions 


n.k  «  0 


(axial  symmetry) 

(no  external  sources) 


(i-1,2,3) 


In  the  direction  of  integration,  these  equations  form  a 
sysLem  of  2.K  equations  for  the  2K  radiances  defined  over 
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a  K-point  linear  mesh.  An  additional  2K  equations  of  similar 
form  are  required  for  each  perpendicular  direction  in  order  to 
define  the  source  function  at  each  grid  point  k.  For  each  path, 
the  transport  equation  yields  a  one-dimensional  solution  matrix 
which  is  a  diagonal  sub-block  of  the  full  three-dimensional 
solution  matrix.  The  individual  paths  are  coupled  through  the 
scalar  source  function  at  intersecting  points  such  that  the  com¬ 
plete  ensemble  of  paths  is  required  to  obtain  closure.  The  complete 
matrix  is  solved  by  iteration  using  the  single-scattering  solution 
(multiple-scattering  sotirce  function  equal  to  zero)  as  the  initial 
guess.  For  the  cases  investigated  to  date  (see  Section  2.6),  two 
to  eight  iterations  (depending  on  the  product  wt)  were  sufficient 

to  achieve  source  function  convergence,  within  two  percent. 

» 

The  importance  of  the  isotropic  source  function  ^transformation 
can  be  appreciated  by  reference  to  the  tangential  (0)  integration 
paths  in  Figure  2-1.  Except  for  the  s-axis ,  the  intersections 
of  these  paths  with  the  cylindrical  shells  are  observed  to  occur  at 
oblique  angles  which  do  not  coincide  with  the  six- angle  quadrature 
of  the  phase  function.  This  situation  presents  no  difficulty  when 
the  source  function  is  isotropic.  However,  a  non-uniform  angular 
quadrature  would  be  required  to  represent  non-isotropic  scattering 
into  the  line-of- sight .  Such  an  approach  is  somewhat  more  complex 
than  the  present  method,  but  offers  the  potential  for  greater 
accuracy  through  higher  resolution.  This  capability  is  provided 
by  the  multiple-flux  model  described  in  the  following  section. 
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2 . 5  Multiple-Flux  Exact  Approach 

The  basis  for  the  exact  approach  developed  under  Task  IA 
is  the  approximation  of  the  multiple-scattering  source  function 
by  an  angular  quadrature  of  arbitrary  resolution.  The  accuracy 
of  the  predicted  radiation  field  increases  as  the  resolution  of 
the  angular  quadrature  increases.  Exact  solutions  are  obtained 
by  increasing  the  resolution  in  successive  calculations  until  con¬ 
vergence  is  achieved. 

The  arbitrary  angular  resolution  of  the  multiple-flux  (hence¬ 
forth  called  "N-flux")  model  affords  considerable  flexibility  in 
the  analysis  of  multiple-scattering  problems.  The  resolution  in 
the  elevation  and  azimuth  planes  is  independently  specified  so 
that  individual  optima  can  be  established.  An  important  feature 
of  the  formulation  is  the  introduction  of  non-uniform  azimuthal 
resolution  which  increases  as  the  radius  increases  in  order  to 
achieve  uniform  accuracy  from  the  isotropic  center  to>  the  aniso¬ 
tropic  edge.  An  alternative  formulation  based  on  uniform  azi¬ 
muthal  resolution  was  considered  early  in  the  model  development 
effort,  but  was  abandoned  due  to  its  greater  complexity  and  lesser 
accuracy  compared  to  the  present  approach. 

Figure  2-4  shows  end-view  and  side-view  projections  of 
representative  lines-of-sight  through  an  assembly  of  infinite 
cylindrical  shells  in  which  the  source  function  for  multiple 
scattering  is  assumed  constant.  In  the  present  development. 
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the  source  function  for  emission  is  also  assumed  constant  over 
the  same  region,  but  this  is  not  a  general  requirement.  The 
unknowns  are  the  radiances  at  the  intersections  of  the  cylindrical 
shells  with  the  parallel  lines-of-sight.  The  constraints  are  the 
transport  equations  for  the  intervening  paths  together  with  the 
boundary  conditions  for  each  line  of  sight.  Together  they  consti¬ 
tute  a  determinate  linear  set.  In  the  figure,  the  rectangular 
z,  y,  x  coordinates  denote  the  plume  axis,  the  line-of-sight 
direction,  and  the  liiie- of- sight  spacing,  respectively.  The 
spherical  r,  0,  <J>  coordinates  denote  the  radius,  azimuth,  and 
elevation  of  the  local  radiance  vectors.  The  integration  paths 
are  determined  by  prescribed  values  of  x  and  <f>  where  the  x's  are 
chosen  to  correspond  to  the  radii  of  the  cylindrical  shells. 

The  integrations  are  performed  over  half  of  the  azimuth  plane 
(-tt/2  <Q<  tt/2)  and  half  the  elevation  plane  (0  <4><  tt/2)  or  one- 
fourth  of  the  radiating  volume  (  tt  steradians)  .  The  remaining 
volume  is  described  by  symmetry  relations  through  simple 
reflections  in  0  and  <p . 

Figure  2-5  shows  the  relationship  between  the  radiances  and 
the  source  functions  at  different  radial  stations  which  provides 
the  coupling  mechanism  for  multiple  lines-of-sight.  This  relation¬ 
ship  is  based  on  the  assumption  of  infinite  axial  symmetry  such 
that  the  properties  of  the  medium  vary  only  with  radius.  Positive 
or  negative  radiance  vectors  are  constructed  from  line-of-sight  rad¬ 
iances  directed  outward  or  inward,  respectively.  In  the  elevation 
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a)  Elevation  Plane 


Figure  2-5a.  Radiances  and  Source  Functions 


gure  2 -5b.  Radiances  and  Source  Functions 


plane.  l  he  sonri'p  function  is  approvi  mater]  bv  an  angular  quadra  - 
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angle.  In  r.he  azimuth  plane,  the  source  function  is  approximated 
by  an  angular  quadrature  of  variable  resolution  which  increases 
in  proportion  to  the  prescribed  number  of  cylindrical  shells. 

The  resolution  increases  from  the  center  (where  the  radiation 
field  is  isotropic  by  symmetry)  to  the  edge  (where  the  radiation 
field  is  most  anisotropic)  such  that  the  accuracy  of  the  quadrature 
is  maintained  independent  of  radius. 

Figure  2-6  shows  the  indexing  system  for  the  matrix  elements 
in  the  azimuth  plane.  The  indices  i  and  k  designate  the  spacing 
(x)  of  the  line-of-sight  and  the  radius  (r)  of  the  cylindrical 
shell.  The  index  j  (which  does  not  affect  the  azimuth  plane 
projection)  designates  the  elevation  angle  ( <p )  of  the  line-of- 
sight  relative  to  the  cylinder  axis.  The  i,k  matrix  elements 
form  an  asymmetric  square  array  in  which  the  ordering  of  the 
subscripts  indicates  the  direction  of  the  path  (outward  or 
inward) . 

Figure  2-7  shows  the  geometric  relations  for  each  matrix  ele¬ 
ment  in  the  azimuth  plane.  The  inclination  of  the  radiance  vec¬ 
tor  relative  to  the  radial  direction  is  given  by  the  complement 
of  the  angle  6ik  or  0ki  where 
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ik 
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ki 


15  cos 
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lk 


=  TT  /  2 


The  line-of-sight  distance 
point  (x-axis)  to  the  point,  k  is  given  by 


or  sk j  ^  measured  from  the  rnid- 
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ndexing  in  Azimuth  Plane 
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For  uniform  properties  (p^r  *,  )  between  points  k  and  k  +  1,  rhe 
transmission  is  given  by 

Tijk  -  exp[-pkKk(sijk+i  -  Eijk )]  • 

With  these  definitions,  the  equation  of  transport  between 
points  k  and  k  +  1  in  the  lower  quadrant  and  points  k  -  1  and  k 
in  the  upper  quadrant  may  be  xrritten  as 
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kji  k+lji  ljk 


to. 
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source  function  8  accounts  tor  multiple  scattering  from  ail 
directions  into  the  line  of  sight. 
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where 


k  =  number  of  discrete  azimuth  angles  at  radius  r. 
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J  “  total  number  of  discrete  elevation  angles 
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phase  function  for  scattering  from  direction 
8i’k’  ‘  anc*  ^t;s  counterparts  into 


direction  0^,  d> j 


solid  angle  for  radiance  vector  or 


subscript,  for  axial  component  (<}> 


0) 


The  axial  component  of  the  source  function  (subscript  zero) 
is  excluded  from  the  double  sum  for  reasons  of  symmetry.  Symme- 
trv  with  respect  to  z  yields  the  following  expression  for  axial 
radiance 


Nok  -  (l  -  “k)Nk°  +  “kSck  • 

Symmetry  with  respect  to  6  yields  the  following  expression  for 
axial  source  function 


J  k 

S  ,  =  Y^P.  7  /n,  +  N.  .  . .  )ft.  +  P  N  ,  ft 

ok  JLmJ  joiL^y  ijk  kji/  xjk  oo  ok  o 

j-1  i“l 

Substitution  of  SQ,C  into  Nq^  and  into  and  elimin¬ 

ates  N  ■  from  the  radiance  matrix.  These  values  are  computed 
independently  subsequent  to  solution  of  the  matrix. 

The  solid  angles  are  evaluated  in  terms  of  the  discrete 
spherical  coordinates  0ik  ^  shown  in  Figure  2-8.  The  solid  angle 
boundaries  are  assumed  to  bisect  the  angular  region  between 
adjacent  radiance  vectors  such  that,  for  an  interior  element, 
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Bor  an  exterior  element,  the  solid  angle  is  one-half  (edge)  or 
one-fourth  (corner)  the  value  given  by  the  above  formula;  i.e., 


i  =  1  j  4  J  (edge) 


fikjk 


-( 


ekk  ~  9k-l,k 


i  =  k  ;  j  /  J  (edge) 


f<iJk 


.  co,(1j  +  *j-i' 


J  - 1 1 

7. — r  cos 

j  =  J  ;  i  /  1 , k  (edge) 


filJk 


-  -(!l*^-)[co3(^l.)  -  COS  ♦/ 


i  =  1  ;  j  =  J  (corner) 


n 


kJk 


-  -(^  -P^-)  COs(^^i)  -  COS 


i  =  k  ;  j  =  J  (corner) 

For  the  axial  component,  th"  solid  angle  ic  given  by 


“  ir(l  -  cos  . 

The  phase  function  is  evaluated  in  terms  of  the  angle  6 
formed  between  the  given  radiance  vector  and  the  line-of -sight . 
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This  angle  is  readily  evaluated  from  the  vector  dot  product  of 
the  corresponding  unit,  vectors  and  r.^^  where 

^ijk  “  *  si-tl  cos  0ip  +  y  uin  j  sin  0ik  +  z  cos  <pj 

Thus , 

Vj'k  '  hjk  '  sln  V  sin  »j[cos(Vk  -  °ik)] 

+  cos  tp .  ,  cos  rf> . 

J  J 

-  cos  8i'j'kij  • 

The  phase  function  is  defined  to  include  scattering  not  only 
from  the  direction  *fc, ‘J’j  »  but  also  from  the  symmetric  direc¬ 
tions  t r  -  0.,.  and  tt  -  <p .  ,  such  that 
1  K  J 

hi'kjj-  ”  p[e(efk-  *j-  Btk’  *j)] 

+  p[8p  -  0i.k,  ♦j..  0ik,  ♦j)] 

+  p[e(9i.k.  »  -  *y.  eik.  Kj)] 

+  p[e(.  -  et,k,  ir  -  +J . ,  eik.  *j)]  ■ 


For  the  axial  component,  the  phase  function  is  given  by 

Pjo  -  2[P(*j)  +pt  -  ♦j)]  ' 

In  this  manner,  the  computed  radiation  field  for  ir-steradians  is 
made  to  account  for  scattering  from  4ir-steradians . 

The  preceding  equations  may  be  solved  by  either  iterative 
or  inversion  methods.  In  the  iterative  approach,  the  equations 
are  suitable  for  programming  without  further  manipulation.  At 


each  iteration,  the  radiances.  Nok,  and  are  computed 

based  on  assumed  values  for  S^,  and  from  the  preced¬ 
ing  iteration  subject  to  the  boundary  condition  »  0  for 

all  i,j  (no  external  sources);  the  assumed  values  for  the  first; 
iteration  are  taken  equal  to  zero  (i.e.,  single  scattering).  In 
the.  inversion  approach,  the  equations  are  manipulated  to  elimin¬ 
ate  the  source  functions  S  in  favor  of  the  radiances  N  which  are 
subsequently  determined  by  matrix  inversion. 

The.  manipulations  required  for  the  inversion  method  arc 
summarized  below.  Elimination  of  in  favor  of  Nqjj.  yields 


Nok(l  “  Wl-? _ !  =  (1  -  oi. 


Ul  A-  Ob  i 

k  oo  o) 


k)Nk°  +  “kL  Epio(»ijk 

i- 1  j  =  i  ' 


-I-  N, 


'kji)fiijk 

Substitution  of  the  Nc^  into  yields 
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Elimination  of  S.^  and  i  in  favor  bf  N^jk  and  ^  yields  a 

KJ (K  -  1)  dimensional  linear  equation  set,  with  KJ (K  -  l)/2 
equations  for  each  quadrant. 


Positive  Quadrant: 
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Negative  Quadrant: 
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(i-l.k  ) 
j  =  1-J 

(k  =  1 ,K-l) 


where. 


*ijk  3  t(x  -  Tijk) 

“ti'kjj'  5  (pii'kjj'  +  xjk“kPjo)!!l,j'k  • 

The  above  equations  are  subject  to  KJ  boundary  conditions 
Nrj1  -  0  (i  -  1,K  ;  j  -  1,J} 

The  above  relations  suffice  to  determine  the  entire  K2J  dimen¬ 
sional  set  of  radiances  for  the  cylinder.  Repeated  application 
at  successive  axial  stations  yields  the  complete  radiation  field 
for  the  plume. 

2 . 6  Computed  Results 

Extensive  parametric  calculations  were  undertaken 
in  order  to  validate  and  compare  the  six- flux  and  N-flux 
scattering  approaches .  The  calculations  included  comparisons 
with  the  predictions  of  other  codes:  the  NASA  finite-difference 
code  for  plane-parallel  geometries  and  the  GAC  Monte-Carlo 
code  for  arbitrary  geometries.  The  calculations  were  per¬ 
formed  for  specified  values  of  the  monochromatic  scattering 
parameters  (optical  depth,  scattering  albedo,  phase  function) 
rather  than  for  a  specified  spectral  region  or  particle  size/ 
density  distribution.  This  procedure  enabled  the  scattering 
codes  to  be  exercised  over  the  full  range  of  optical  environ¬ 
ments  without  the  attendant  complexities  of  a  complete  system 
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simulation.  The  evaluation  of  the  optical  parameters  In  terms 
of  the  physical  properties  of  a  specific  gas/partiole  medium  is 
considered  elsewhere  in  this  report. 

2.6.1  Planar  media,  The  relative  simplicity  and  general 
availability  of  solutions  for  planar  media  make  such  solutions 
ideal  for  the  preliminary  test  of  physical  models  and  parametric 
effects.  This  section  summarizes  results  obtained  for  uniform 
one-  and  two- dimensional  planar  media  using  an  explicit  six-flux 
model  (Method  I;  Figures  2-9  through  2-12)  and  a  standard  finite- 
element  program  (NASA  code;  Figures  2-13  through  2-15)  modified 
to  include  thermal  emission. 

Figure  2-9  shows  the  transmissivity,  reflectivity,  and  emis- 
sivity  of  a  unifora  planar  medium  of  unit  optical  depth  at  normal 
aspect  for  various  phase  functions  (isotropic,  fore  and  aft, 
forward  only)  In  addition,  isotropic  scattering  solutions  are 
compared  for  one-  and  two- dimensional  planar  media  (slab  vs. 
ribbon)  in  order  to  bound  the  possible  range  of  values  for 
cylindrical  media.  In  all  cases,  scattering  is  observed  to 
reduce  the  emitted  component  and  to  increase  the  transmitted- 
plus-ref lected  component  of  the  radiation  signature.  Comparisons 
among  the  slab  solution  for  different  phase  functions  show  the 
effects  of  backscattering  (b>0)  and  side-scattering  (s>0)  on 
the  reflected  and  emitted  components.  Comparisons  between  the 


2-51 


uirixpri 


nnwwnw* . .  ii«  "»  ^  •uremiimi 


slab  and  ribbon  solutions  for  isotropic  scattering  (f“bMs) 
reveal  essential  differences  between  conservative  (e  +  t  +  P  »1) 
and  non- conservative  (c  +  i  4-  p  <1)  radiation  fields  which  would 
not  be  predicted  by  a  two-flux  model  (s«0). 

Figure  2-10  shows  the  transmissivity,  reflectivity,  and  emis- 
sivity  of  a  uniform  slab  at  tangent  aspect  for  the  same  range 
of  parameters  as  in  Figure  2-9.  In  this  case,  the  incident  ra¬ 
diation  is  assumed  to  be  perpendicular  to  the  surface  while  the 
observed  radiation  is  parallel.  Thus,  the  transmitted  radiation 
component  undergoes  a  change  in  direction  from  the  incident  side 
to  the  emergent  side  which  accounts  for  t  0  as  co  -+■  0,  The 
emitted  component  parallel  to  the  slab  is  black  (infinite  optical 
depth)  cnly  in  the  limit  co  0;  for  finite  m,  the  emissivity  is 
reduced  due  to  losses  through  the  slab  boundaries.  Note  that 
this  coupling  between  the  normal  and  tangential  radiance  components 
does  not  occur  for  the  tvc-fiox  models. 

Figure  2-11  scows  the  radiance  of  a  uniform  emitting  planar 
medium  at  normal  aspect  over  a  range  of  optical  depths  in  the 
absence  of  external  radiation  sources,  The  values  are  normalized 
by  the  non-scattering  solution  for  co  =  0  and  are  therefore 
equivalent  to  normalized  emissivities .  The  range  of  scattering 
parameters  and  geometries  is  the  same  as  in  Figures  2-9  and 
2- 10.  In  all  cases,  increased  scattering  (<o)  is  observed  to 
reduce  the  observed  radiance  for  constant  x  due  to  equal  but 
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Figure  2-11.  Normal  Radiance  versus  Albedo 
for  Uniform  Planar  Medium 
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opposite  variations  in  the  optical  depths  for  absorption  and 
scattering.  The  effect  of  side- scattering  is  observed  to 
increase  the  radiance  from  the  slab  relative  to  the  ribbon 
and  from  the  six-flux  relative  to  the  two-flux  (f=b  in  both  cases) 
but  to  decrease  the  radiance  at  large  optical  depths  relative 
to  the  forward-scattering  case  (f“l). 

Although  not  shown,  the  normalized  axial  radiance  component 

(Nz  (oj) /Nz  (0)]  exhibited  close  similarity  to  the  normalized 

radial  component  [n^.  (to) /N^  (0)]  for  isotropic  scattering  but 

maintained  constant  value  [N  (co)/N  (0)  =  l]  for  two-stream 

z  z 

scattering.  The  normalizing  values  N^.(O)  and  N  (0)  are  simply 

the  given  black  body  function  (N°)  times  the  apparent  emissivity 
-T 

(1-e  or  unity,  respectively).  Thus,  the  six-stream  approxi¬ 
mation  accounts  for  the  redistribution  of  radiance  (from  axial 
to  radial)  due  to  scattering  whereas  the  two- stream  approxi¬ 
mation  does  not. 

Figure  2-12  shows  analogous  results  to  Figure  2-11  except 
that  the  optical  depth  parameter  is  based  on  the  absorption 
component  rather  than  the  total  extinction  (absorption  4  scattei 
ing) .  In  this  case,  the  emitted  radiation  remains  constant  as 
a)  increases  rather  than  decreasing,  as  before,  in  proportion  to 

(1-ui).  As  a  result,  the  curves  lie  above  those  in  Figure  2-11 
by  an  amount  which  increases  with  increasing  w.  In  the  case  of 
forward  scattering,  the  radiance  remains  finite  in  the  limit 

u)  -+■  1 . 
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Figure  2-13  shows  a  comparison  between  the  predictions  of 
the  two- flux  and  six- flux  models  and  the  NASA  finite-element 
code  for  the  normalized  radiance  at  normal  aspect  from  an  iso- 
tropically-scattering  uniform  slab  of  specified  optical  depth. 

For  small  optical  depth  (t**.1),  the  NASA  code  is  in  close  agree¬ 
ment  with,  the  two-flux  model  (f=b“^).  For  large  optical  depth 
(ri!10>  ,  the  NASA  code  is  in  closo  agreement  with  the  six-flux 
model  (f~b=s=l/G) .  For  intermediate  optical  depth  (based  on  a 
single  NASA  code  calculation) ,  the  NASA  code  prediction  is  inter¬ 
mediate  between  the  other  two.  Taking  the  NASA  code  as  the 
standard,  this  comparison  provides  physical  insight  into  the 
sources  of  error  which  can  result  from  coarse  angular  quadrature. 
In  short,  the  twe-  and  six- flux  quadratures,  respectively,  under - 
predict  and  over-predict  the  effective  emissiv.i.ty  of  the  medium. 
However,  from  the  comparative  results  of  Figures  2-9  to  2-13  for 
slabs  and  ribbons,  the  error  in  the  six-flux  approximation  at 
small  optical  depths  should  be  appreciably  reduced  for  cylindrical 
geometries  (i.e.  plumes  ). 

Figures  2-14  and  2-15  show  further  comparisons  between 
•  predictions  of  the  six-flux  model  and  the  NaSA  code  for 
an  isotropically  scattering  uniform  slab.  These  figures  show 
the  variation  with  aspect  angle  of  the  radiance  and  radiant 
intensity  (radiance  times  projected  area)  normalized  to  their 
respective  values  at  normal  aspect.  The  predictions  of  the  six- 
flux  model  (Method  I)  represent  interpolations  between  the  values 
actually  computed  for  normal  and  tangent  aspect  and,  therefore. 
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have  limited  quantitative  accuracy  between  those  limits . 

However,  the  predictions  and  comparisons  are  useful  to  illustrate 
trends.  In  particular,  both  models  predict  che  same  qualitative 
variation  of  radiance  with  aspect  (Figure  2-14)  which  is  strongly 
dependent  on  optical  depth.  For  small  optical  depths  and  low 
scattering  albedos,  the  radiance  increases  toward  tangent 
aspect  due  to  the  increased  length  of  the  optical  path.  However, 
for  large  optical  depths  and  high  scattering  albedos,  the  radiance 
decreases  toward  tangent  aspect  due  to  side-scattering  losses 
not  predicted  by  two- flux  models.  The  aspect  dependence  of  the 
radiant;  intensity  (Figure  2-13)  includes  the  variation  in  projected 
area  (proportional  to  sin  <f>)  and,  consequently,  shows  a  monotonic 
decrease  from  normal  to  tangent  aspect.  In  both  Figures,  the 
six- flux  model  predicts  lower  levels  of  radiation  for  higher 
scattering  albedos,  with  greater  sensitivity  for  larger  optical 
depths.  The  opposite  behavior  of  the  NASA  predictions  for  low 
optical  depth  does  not  lead  to  larger  absolute  differences. 

2.6.2  Cylindrical  media.  Compared  to  the  vast  literature 
on  multiple-scattering  solutions  for  planar  media,  relatively 
little  work  has  been  published  for  cylindrical  media.  Accounting 
for  this  essential  geometric  effect  was  a  central  feature  of  the 
SIRRM  model  development.  Hence,  the  logical  first  test  of  a  new 
solution  for  cylindrical  media  is  a  comparison  with  previous 
solutions  for  planar  media. 
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Figure  2-16  shows  the  normal  radiance  of  a  uniform  infinite 
cylinder  (centered  line-of-sight)  and  of  a  uniform  slab  of  equal 
optical  depth.  The  predictions  are  based  on  the  six-flux  model 
(Method  I  or  II)  and  on  the  GAC  Monte  Carlo  code.  Except  in  the 
non-scattering  limit  (gj=0)  ,  the  radiances  are  observed  to  be 
lower  for  the  cylinder  than  for  the  slab  due  to  side-scattering 
losses  (not  predicted  by  two-flux  models) .  The  six-flux  and  Monte 
Carlo  results  are  shown  to  be  in  essential  agreement  for  both 
isotropic  and  forward  scattering  phase  functions.  The  forward 
scattering  result  is  only  slightly  greater  than  the  isotropic 
scattering  result  due  to  the  assumption  of  uniform  properties; 
larger  sensitivities  would  be  expected  for  the  non-uniform  con¬ 
ditions  of  actual  plumes. 

Figure  2-17  shows  the  normal  radiance  for  isotropic  scattering 
of  a  uniform  infinite  cylinder  (centered  line-of-sight  )  and  of 
a  uniform  slab  and  ribbon  of  equal  optical  depth.  The  predictions 
are  based  on  the  6-flux  model  (Method  I  or  II)  with  comparative 
calculations  from  the  N-flux  model  and  the  Monte  Carlo  code. 

The  variation  of  radiance  with  optical  depth  (upper  left  plot) 
is  qualitatively  consistent  with  previous  predictions  for 
planar  media  (Figure  2-11).  Specific  comparisons  for  equal 
optical  depths  (remaining  plots)  show  a  uniform  progression  of 
the  cylinder  radiances  from  the  ribbon  values  at  small  optical 
depths  to  the  slab  values  at  large  optical  depths.  This  be¬ 
havior  is  the  result  of  side-scattering  losses  from  the  cylinder 
which  increase  with  decre: sing  optical  depth.  These  losses 


Figure  2-17.  Normal  Radianc 


result  in  improved  predictive  accuracy  for  the  six-flux  model 
it  small  optical  depths.  Comparison  of  the  six-flux  model  pre¬ 
dictions  with  the  more  exact  predictions  of  the  N-flux  model  and 
Monte  Carlo  codes  reveals  a.  residual  over-prediction  for  small 
optical  depths  but  substantial  improvement  relative  to  similar 
predictions  for  a  slab  (compare  Figure  2-13')  .  Comparison  of  the 
N-flux  and  Monte  Carlo  predictions  for  unit  optical  depth  revels 
excellent  agreement.  The  comparison  is  close,  but  less  exact  at 
an  optical  depth  of  ten. 

Figure  2-18  shows  the  tangent  radiances  (i.e.,  aspect  =  0 
or  180  degree?)  for  the  same  conditions  as  the  normal  radiances 
in  Figure  2-17.  As  before,  the  curves  are  shown  normalized  by 
their  values  for  a  non- scattering  medium,  i.e.,  the  local  Planck 
function  [  Nz  (oj=0)=N°]  .  The  curves  reveal  the  same  qualitative 
trends  as  the  normal  radiances  in  Figure  2-17  except  for  a  reduced 
sensitivity  to  optical  depth  due  to  cylindrical  curvature.  As  n 
result,  the  tangent  radiances  remain  close  to  the  ribbon  values 
for  all  optical  depths.  This  effect  is  confirmed  by  the  more 
exact  predictions  of  the  N-flux  model  and  Monte  Carlo  code. 

Figure  2-19  shows  the  aspect  dependence  of:  the  eirassi  vity 
of  a  uniform  infinite  cylinder  (centered  line-of-sight)  as 
predicted  for  a  given  albedo  by  two- flux  and  six-flux  model 
approximations  (Method  II) .  The  curves  are  qualitatively  consis¬ 
tent  (but  quantitatively  different)  for  small  optical  depths  near 
normal  aspect.  However,  for  large  optical  depths  near  tangent 
aspect,  the  two  approximations  prodirt  opposite  trends  due  to 
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the  effect  of  wide-scattering  losses  ne&r  grazing  incidence. 

This  behavior  is  consistent  with  previous  predictions  of  radiance 
vs.  aspect  for  a  planar  slab  (Figure  2-14)  and  may  have  signifi¬ 
cant  negative  impact  on  the  detectability  of  plumes  from  nose 
aspect.  Comparison  of  the  predictions  of  Method  I  and  II  for 
normal  and  tangent  aspect  reveals  good  agreement  for  both  phase 
function  approximations  at  comparable  optical  depths.  Comparison 

of  the  predictions  of  Method  II  ano  the  Monte  Carlo  code  for 
isotropic  scattering  at  unit  optical  depth  reveals  agreement  within 

20  percent  for  all  aspect  angles. 

Figure  2-20  shows  a  comparison  between  the  six-flux 
(Method  I)  and  N-flux  models  for  the  normalized  radiance 
from  an  isotropically  scattering  uniform  infinite  cylinder 
(centered  line-of- sight)  over  the  full  range  of  albedo  and 
optical  deoth.  These  curves  are  equivalent  to  normalized 
emissivities  from  Figure  2-19.  As  in  Figure  2-3.4  (slab),  the 
present  six- flux  predictions  represent,  interpolations  between 
the  values  actually  computed  for  normal  and  tangent  aspect  and, 
therefore,  have  limited  quantitative  accuracy  between  those 
limits.  In  particular,  the  Method  I  interpolation  formulae  are 
observed  to  under-predict  relative  to  the  N-f lux  model  but  to 
yield  valid  qualitative  trends.  Both  models  yield  comparable 
predictions  at  normal  and  tangent  aspects.  Comparison  of  N-flux 
and  Monte  Carlo  predictions  are  in  excellent  agreement  and 
appear  to  confirm  the  6-flux  Method  II  predictions  in  Figure 
2-19. 
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Figure  2-21  shows  the  N-flux  model  prediction  for  the 
variation  of  normal  radiance  with  distance  from  the  cylinder 
axis  for  the  same  conditions  as  Figure  2-20.  Integration  under 
these  curves  would  yield  the  normalized  station  radiation  at 
side  aspect.  The  effect  of  scattering  is  to  reduce  the  radiance 
from  the  off-axis  regions,  with  greater  sensitivity  for  larger 
optical  depths.  The  assumption  of  the  six-flux  model  I  that  the 
radiance  profile  is  independent  of  albedo  is,  therefore,  seen 
to  be  in  error  for  large  optical  depth.  This  assumption  is  not 
made  by  six-flux  model  II. 

Figure  2-22  shows  a  comparison  between  the  six-flux 
(Method  I)  and  N-flux  model  predictions  for  station  radiation 
as  a  function  of  aspect  angle  for  the  same  conditions  as  Figure 
2-20.  Both  models  predict  the  same  qualitative  effect  of 
increasing  albedo,  namely,  decreased  radiation  relative  to  side 
aspect  with  greater  sensitivity  fcr  large  optical  depth.  This 
behavior  is  consistent  with  previous  6-flux  model  predictions 
for  a  planar  slab  (compare  Figure  2-15)  except  that  the  sensitivity 
to  albedo  is  greater  in  the  present  case.  However,  the  six- 
flux  model  (Method  I)  interpolation  formulae  appear  to  over¬ 
predict  the  sensitivity  as  predicted  by  the  N-flux  model  and  as 
confirmed  by  additional  predictions  with  the  GAG  Monte  Carlo 
code.  The  latter  two  predictions  are  shown  to  be  in  excellent 
agreement . 
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The  preceding  results  for  infinite  cylindrical  media 
represent  valid  solutions  for  rocket  plumes  provided  that 
the  variation  in  axial  radiance  is  negligible  over  a  distance 
of  one  optical  depth  (i.e.,  dN  /dT  <<N  ).  This  criterion  is 
satisfied  throughout  most  of  the.  plume  mixing  region  but  may  not 
be  satisfied  in  local  regions  of  the  plume  inviscid  core  when 
strong  shocks  are  present.  In  order  to  assess  the  impact  of 
shock- induced  gradients  on  the  coupled  gas/particle  radiation  field, 
the  six-flux  model  (Method  II)  was  extended  to  finite  cylindrical 
media  with  arbitrary  end-wall  boundary  conditions.  Numerical 
experiments  were  recently  initiated  to  assess  the  impact  of  a 
.  i  ong  extra*  nai  ioiii.at.ion  source  parallel  to  t  he  axis  on  the 
radiation  scattered  to  side  aspect.  A  key  objective  of  this 
assessment  Is  the  determination  of  the  axial  region  of  influence 
beyond  which  the  infinle  cylinder  model  becomes  a  valid  approx 

n. "  i i mi 


Figure  2-23  shows  the  result  of  a  preliminary  test  cal¬ 
culation  with  the  extended  six-flux  model.  In  this  case,  the 
'  i  n  ,nr  is  assumed  uni  form  with  no  external  radiation 
■  •  •  rt  the  computed  radiation  field  is  symmeti  it  about  . 

cylinder  mid-point  The  radial  and  axial  emissivities  are  observed 
to  decrease  from  the  mid-points  toward  the  corners  as  a  result 
of  three-dimensional  radiation-coupling  effects.  For  the  giver, 
optical  depths,  the  mid-point  radial  emissivity  is  nearly  equal 
to  the  emissivity  of  an  infinite  cvlmder;  i.e.,  the  mid-point 
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values  are  nearly  independent  of  end-wall  boundary  conditions . 
Thus,  the  axial  region  of  influence  is  contained  within  the  finite 
cylinder  considered.  Comparative  Monte  Carlo  predictions  appear 
to  validate  the  six-flux  model  result. 

As  a  last  test  case,  a  condition  of  practical  significance 
was  chosen.  It  involves  an  occulted  emitting  core  and  a  non- 
emitting  scattering  envelope,  simulating  a  nose  aspect  signature 
(see  Figure  2-24).  This  is  a  stringent  test  case  because  the 
nose  aspect  signature  is  due  entirely  to  scattering  and  because 
the  core  is  a  localized  source  leading  to  a  highly  anisotropic 
rr.dia:  ion  field.  A  comparison  between  the  forward  scattered 
radiance  predicted  by  the.  6-flux  model  (Method  II)  and  the  Monte 
Carlo  method  showed  agreement  within  20%. 


Figure  2-24.  Core/Envelope  Structure  for 
Nose  Aspect  Test  Case. 
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Summary  of  Models  and  Results 


The  results  discussed  in  Section  2.6  include  calculations 
for  planar  and  cylindrical  media  using  6-flux,  finite-difference, 
and  Monte  Carlo  solution  techniques.  The  calculations  include 
the  full  range  of  optical  parameters  [  albedo  (0  to  1) ,  optical 
depth  (.1  to  10),  phase  function  '(isotropic  to  delta  function)] 
and  observables  (transmissivity,  reflectivity,  emissivity; 
radiance  and  station  radiation  versus  position  and  aspect) .  The 
results  demonstrate  the  importance  of  three-dimensional  multiple¬ 
scattering  effects  and  the  relative  ability  of  different 
physical  models  to  provide  an  accurate  simulation. 

Three  new  models  were  developed  and  verified  in  this  phase 
of  the  effort.  The  first  was  an  explicit  six-flux  model  which 
is  exact  (for  six- flux,  phase  functions)  for  planar  media  and 
approximate  for  cylindrical  media.  The  second  was  an  implicit 
six-flux  model  which  is  exact  (for  six-flux  phase  functions)  for 
cylindrical  media,  The  third  was  an  implicit  multiple -flux 
("N-flux")  model  which  is  exact  (for  arbitrary  phase  functions) 
for  cylindrical  media.  The  six-flux  models  provide  the  minimum 
angular  resolution  necessary  to  account  for  three-dimensional 
scattering  effects.  As  such,  they  constitute  valid  engineering 
approximations  which  are  highly  efficient  for  production  calcu¬ 
lations.  The  N-flux  model  provides  complete  flexibility  in  the 
selection  of  angular  resolution  (elevation  and  azimuth)  .  As  r. 
result,  alternate  levels  of  resolution  can  be  tested  and  compared 
with  respect  to  accuracy  and  efficiency. 
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The  six-flux  models  are  distinguished  principally  by  the 
degree  of  rigor  in  treating  side-scattering  effects.  That  is, 
Method  I  approximates  these  effects  via  simple  closure  relations 
whereas  Method  II  computes  these  effects  via  exact  coupling 
equations.  Both  methods  solve  initially  for  the  normal  and 
tangential  radiances  and  subsequently  for  the  radiance  at 
arbitrary  aspect  (if  required) .  The  fundamental  difference  is 
associated  with  the  initial  solution  step. 

To  the  authors'  knowledge,  the  N-flux  model  represents 
the  first  application  of  general  finite-element  techniques  to 
emitting  cylindrical  media.  Previous  application  of  this  tech¬ 
nique  was  to  non-emitting  planar  media  (Barkstrom,  1976).  The 
information  content  of  the  calculation  is  considerable.  For 
each  radial  station,  the  model  predicts  the  complete  three- 
dimensional  radiance  distribxition  in  spherical  coordinates  . 

Typical  results  for  two  uniform  isotropically  scattering  cy¬ 
linders  of  different  optical  depth  are  illustrated  in  Figure  2-25. 
The  figure  shows  polar  plots  in  azimuth  and  elevation  of  the  local 
radiance  distribution  at  a  series  of  radial  stations  from  the 
c  ?>  ter  of  each  cylinder  to  the;  edge.  In  each  case,  the  rat':’  : 
are  normalized  by  the  edge  value  for  a  centered  line-of-sight 
at  side  aspect  (0*<}>=Tr/2)  .  In  the  azimuth  plane,  the  radiance 
distribution  is  observed  to  be  isotropic  at  the  center  and  to 
evolve  toward  a  forward  lobe  at  the  edge.  In  the  elevation 
plane,  the  radiance  distribution  is  roughly  eliptic  at  the  center. 
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(approaching  isotropic  for  large  optical  depths)  and  evolves 
toward  a  composite  forward/ side  lobe  at  the  edge.  The  radiances 
are  observed  to  increase  or  decrease  from  the  center  to  the  edge 
for  small  or  large  optical  depths,  respectively,  with  the  maximum 
rate  of  change  occuring  near  the  surface.  This  result  may  be 
compared  to  a  non-scattering  medium  in  which  the  radiance 
always  increases  monotonically  outward. 

The  N-flux  model  provides  a  powerful  tool  for  the  analysis 
of  plume  scattering  problems  of  arbitrary  complexity.  The 
accuracy  of  the  prediction  is  limited  only  by  the  computational 
resources  of  the  user.  Because  of  the  significant  cost  im¬ 
plications  of  repetitive  calculations  for  multiple  spectral 
points,  an  important  aspect  of  the  continuing  effort  is  the 
optimization  of  the  mathematical  procedures  and  the  grid  selection 
(spatial  and  angular).  Mathematical  optimization  includes  the 
selection  of  iterative  or  inversion  procedures  for  small  or  large 
scattering  optical  depth  (oix)  and  the  use  of  a  different  (finer) 
radial  grid  for  evaluation  of  the  emission/extinction  terms. 

Grid  optimization  involves  the  correlation  of  accuracy  and 
efficiency  with  both  the  number  of  elements  (radius  and  aspect) 
and  the  positions  of  their  boundaries. 

A  rough  preliminary  indication  of  the  degree  of  resolution 
required  for  a  given  accuracy  is  shown  in  Figure  2-26.  This 
figure  is  representative  of  the  limited  computational  experience 


acquired  to  date  but  does  not  necessarily  possess  broader  quan¬ 
titative  applicability.  The  results  shown  are  for  an  optically 
thick  uniform  cylinder  with  isotropic  scattering  (t=10,  a)**. 8). 

No  particular  effort  was  made  to  optimize  the  grid  selection  for 
a  given  number  of  points .  The  accuracy  is  observed  to  be  much 
more  sensitive  to  the  number  of  radial  elements  than  to  the 
number  of  aspect  angles . 

By  far  the  largest  portion  of  the  computing  time  is  in  the 
calculation  of  the  source  functions.  As  a  result,  CPU  time  is 
roughly  proportional  to  the  square  of  the  number  of  aspect 
angles  ( cp  -  4> j  •  )  and  to  the  cube  of  the  number  of  radial  segments 
(rk,  Oi ,  0.,  »  )  ,  Thus,  doubling  both  the  number  of  aspect  angles  and 
the  number  of  radial  segments  results  in  a  computing  time  increase 
of  about  a  factor  of  30.  Optimization  of  the  selection  of  these 
numbers  is,  therefore,  of  considerable  importance  for  production 
calculations . 

Figure  2-27  summarizes  the  parametric  and  comparative  cal¬ 
culations  performed  to  date  as  reported  in  Sections  2.6  and  2.7. 
This  summary  illustrates  the  interrelationships  among  the  various 
techniques  and  :he  scope  of  physical  problems  considered.  The 
emphasis  in  this  phase  of  the  study  was  on  model  development 
and  verification.  Additional  work  is  required  in  the  areas  of 
model  optimization  and  application.  In  particular,  final  quan¬ 
titative  assessment  of  the  relative  accuracy  and  efficiency  of 
the  different  calculation  methods  in  different,  physical  situa¬ 
tions  must  necessarily  await  further  model  refinements  and 


2-83 


'•*1  P -  r^-tPif  i;'i  u  J  Vi^r-T*  Ttj  y  --  - 


additional  parametric  predictions.  These  efforts  are  currently 
underway.  However,  based  on  the  parametric  and  comparative 
calculations  performed  to  date,  certain  general  conclusions  can 
be  drawn: 

1)  The  six-flux  angular  quadrature  (Method  I  or  II)  yields 

quantitatively  accurate  predictions  for  normal  and  tangent 
aspect  in  the  cases  of  intermediate- to- large  optical  depth 
or  anisotropic  scattering.  The  standards  for  comparison 
are  the  NASA  finite-element  code  for  slabs  and  the  N-Flux 
and  Monte-Carlo  codes  for  cylinders. 

2>  The  six-flux  angular  quadrature  (Method  I  or  II)  tends  to 
overpredict  the  normal  and  tangent  radiation  for  isotropic 
scattering  at  small  optical  depths.  The  error  is  greatest 
for  slabs  and  least  for  ribbons  due  to  the  corrective  effect 
of  side-scattering  losses.  For  cylindrical  plume  geometries, 
the  error  appears  to  be  acceptable  (less  than  20  percent). 

3)  The  aspect  dependence  predicted  by  the  six- flux  methods  is 
based  on  prior  solution  for  normal  and  tangent  aspect.  For 
intermediate  aspect  angles,  simple  interpolation  (Method  I) 
gave  the  correct  qualitative  behavior,  but  overestimated 
the  sensitivity  to  albedo  by  comparison  to  the  N-flux  model. 
Line-of -sight  integration  (Method  II)  gave  the  correct 
quantitative  behavior  by  comparison  to  limited  Monte  Carlo 
calculations  for  intermediate  optical  depth. 

4)  The  N-Flux  and  Monte  Carlo  calculations  demonstrated  uni¬ 
formly  excellent  agreement  for  all  cases  investigated. 
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BAND  MODEL  APPROACHES 


3.1  Introduction 


The  equation  of  radiative  transfer,  which  is  rigorous  on  a 
monochromatic  basis,  is  only  approximate  when  applied  over  a 
finite  spectral  interval.  The  spectral  averaging  procedures 
required  for  the  calculation  of  radiant  transfer  in  a  multi¬ 
ple  scattering  medium  must  be  considered. 


The  integral  form  of  the  equation  of  transfer  for  a  finite 
spectral  interval  AA  may  be  written 


-  -  r  o  dT  A<-s 
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where  a  spectrally  averaged  quantity  x  is  defined  as 


XA  "  TX  /  xAdX 
AA 


and  where  the  relatively  small  spectral  variations  of  and 
within  AA  are  neglected.  The  integrand  contains  three  spectrally- 
averaged  terms  involving  the  transmittance  derivative  and  its 


3-1 


Product  with  the  albedo  and  radiance.  Of  these  terms ,  only  the 
transmittance  derivative  is  amenable  to  simple  averaging  through 
superposition  of  the  individual  species  averages.  Evaluation 
of  the  product  terms  requires  additional  approximations  which 
are  investigated  in  this  section.  The  goal  of  this  investigation 
is  to  develop  as  comprehensive  a  band  model  theory  of  combined 
scattering  and  gaseous  absorption  as  possible  and  to  determine 
the  functional  dependence  of  the  mean  transmissivity  and  emissivity 
on  the  gaseous  band  model  and  particle  scattering  parameters 
as  a  guide  to  their  most  reasonable  combination  in  the  more 
general  three-dimensional  numerical  solutions  of  the  radiative 
transfer  equation. 


Tti  developing  a  combined  absorption-scattering  band  model 
formalism,  the  approach  has  been  used  which  is  analogous  to 
that  of  the  classical  gaseous  band-model  theory,  namely,  from  a 
monochromatic  solution  of  the  radiative  transfer  equation,  to 
perform  a  spectral  integration  over  a  defined  gaseous  absorption 
band,  making  approximations  as  required  to  obtain  a  closed  form 
solution  (and  studying  the  limitations  on  the  accuracy  of  the 
solution  imposed  by  the  required  approximations) . 


3.2.  Band  Model  Formulation  For  One-Dimensional  Scattering 

A  starting  point  has  been  selected  to  be  the  "one-dimensional" 
(i.e.,  two-flux)  scattering  formalism  of  Vanderbilt  and  Slack 
(1976),  which  has  the  advantage  of  having  been  developed  to  the 
state  in  which  various  closed- form  solutions  exist  which  are 
explicit  in  certain  limiting  casts. 
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It  is  important  to  note  that  the  two-flux  and  six-flux  ap¬ 
proaches  to  the  problem  of  radiation  from  a  plane-parallel  slab 
yield  similar  functional  forms,  differing  only  in  the  definitions 
of  the  effective  backward  and  forward  scattering  cross  sections. 
Thus,  the  so-called  "one-dimensional"  scattering  solution  is 
capable  of  reproducing  the  effects  of  three-dimensional  scattering 
in  this  particular  case. 

It  is  noted  also  that  the  condition  of  completely  anisotropic 
scattering  which  is  sometimes  evoked  as  justification  for  applica¬ 
tion  of  the  one-dimens J onal  formalism  is  a  sufficient  condition, 
but  by  no  means  a  necessary  condition  for  its  applicability. 
Because  of  its  relative  mathematical  simplicity  (in  addition  to 
the  fact  that  the  previously  mentioned  sufficient  condition  for 
its  validity,  i.e.,  the  assumption  of  complete  anisotropy,  is 
very  closely  approached  in  many  real  scattering  problems), 
this  formalism  will  always  be  of  great  value  in  a  large  number 
of  practical  situations. 


We  start  with  one  approximate  solution  of  the  two-flux 
scattering  equation,  which  is  exact  in  a  limiting  case  and  which 
provides  an  iterative  solution  under  more  general  conditions 
(Vanderbilt  and  Slack,  1976): 
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where 


W  m  ^*Naax^ds' 


(3-2) 


1+ 

I  a. 


-  1+ 


2Bu) )' 

ra 


(3-3) 


In  these  equations,  N  is  the  particle  density,  a  the  absorption 

a 

cross  section  and  u  ^  the  Lackscattering  cross  section,  and  a 
is  the  ratio  of  backscattering  to  total  scattering.  For  simplicity, 
the  absorption  and  scattering  cross  sections  are  referenced  to 
a  single  particle  density;  however,  for  mixtures  of  particles  or 
particles  plus  gases,  the  expressions  are  readily  generalized. 

For  the  case  of  constant  x*  an  exact  solution  is  provided 
for  an  otherwise  inhomogeneous  path,  i.e., 


I'(0)  = 
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where 
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The  absorption  optical  depth  is  denoted  by  t  ,  which  may 
represent  the  sum  of  contributions  from  a  number  of  gaseous 
and  particulate  species. 

For  the  isothermal  case,  direct  integration  provides  (on 
division  by  the  Planck  blackbody  function)  the  known  expression 
for  emissivity: 


On  factoring  and  rearranging  into  a  form  more  suitable  for 
our  purposes ,  we  obtain 


Briefly  considering  the  limit  of  small  pco  (x^l)  ,  we 
obtain  a  somewhat  more  transparent  form 


,1  goo  (,  ,  ~Ts  Ta\ 

1+2  r^r  u+e  1 

Because  of  the  relatively  small  variation  of  the  denominator 
from  unity,  the  mean  value  over  a  spectral  interval  can  be 
written  as 
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where 


(3-10) 


represents  the  band-model  curve-of-growth  for  the  gas  plus  par¬ 
ticle  absorption  component  (that  is,  a  known  function).  Equation 
(3-9)  explicitly  shows  the  modifications  resulting  from  the 
simultaneous  presence  of  scattering. 


For  x  significantly  greater  than  unity,  it  still  appears  rea¬ 
sonable  to  consider  the  preponderant  effect  of  spectral  variation 
to  occur  in  the  numerator  (which  varies  from  unity  to  zero) . 

The  percentage  error  introduced  by  this  assumption  will  be  in¬ 
vestigated  in  the  numerical  studies  in  Subtask  II.  We  write 
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where 
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We  designate  the  numerator  of  Eq.(3-il)  as  1-t1  (since  it 
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corresponds  to  the  emissivity  =  1- transmi ssivity  in  the 
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3.2.1  Elsasser  model.  We  examine  the  behavior  of  the  func¬ 
tion  t'  when  the  absorption  results  from  one  gaseous  component 
represented  by  the  Elsasser  model  (for  which  case  the  effects  of 
the  gaseous  spectral  structure  will  be  the  greatest) : 
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sinh  3 
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where 


=  ku  =  Su/d 
(mean  absorption  optical  depth) 
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(fine-structure  parameter) 


and 
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The  properties  of  this  integral  will  be  examined  in  detail  for 
the  limiting  cases  of  small  3  and  large  6. 


For  large  values  of  the  parameter  6  (i,e.,  for  significant 

overlapping  of  the  spectral  lines),  the  function  t'  becomes 
(on  expansion  of  the  terms  in  l-cosx/cosh8) , 


t'  “  exp  (Ta2tanh28+2iaTgt&nh8)  J 
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where  is  the  modified  Bessel  function  of  the  first  kind, 

o 

1  2 

Since  exp(-x)~l-x  for  small  x,  and  I^y)  ~  1+^y  for 
small  y  (and  y  <<  x) ,  the  I  term  provides  a  higher-order 
correction  to  the  exponential  term.  Thus  a  useful  approximation 
is 

t'  =  exp  [- (7a2tanh2B+2xaTstauh  0  |  (3-19) 

Note  that  for  no  scattering  (xg  =  0) ,  we  have 

t  =  1-7  =  7'  =  exp(-?atanhB) ,  (3-20) 

the  classical  Elsasser  model  result. 

For  the  more  interesting  case  of  8  small  (little  over¬ 
lapping  of  spectral  lines) ,  we  have  (using  cosh  8  ~  1) 
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By  making  the  substitution 
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to  transform  to  a  more  standard  form  (e.g.,  that  of  a  Laplace 
transform) ,  the  function 
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becomes 
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This  form  is  somewhat  unwieldy,  and  does  not  appear  to  be 
immediately  integrable.  For  the  limiting  case  a  -  0  (no 
scattering) ,  the  expression  above  reduces  to 


t'  -  |  f  exp  (-  /Ez)  - - — r-  dz 


z(z-l)' 


=  erfc(b)£)  =  erfc  |(x  8/2) 


the  usual  Elsasser  model  result  for  small  8  . 


In  the  other  limiting  case  of  b  *  0  (more  strictly 
speaking,  b/a  -*■  0)  the  integral  for  t’  reduces  to 
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This  expression  may  be  converted  to  a  known  Laplace  transform, 
yielding 


t' 
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where 
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as  before,  and  Kq(x)  is  the  modified  Bessel  function  of  the 
second  kind.  Such  functions  (of  integral  order)  appear  to  be 
the  least  tractable  of  the  Bessel  functions.  However,  since 
Kq(x)  is  only  weakly  (logarithmically)  divergent  at  the  origin, 
series  expressions  may  be  used  and  integrated  over  the  finite 
interval  indicated. 

Because  of  the  mathematical  intractability  of  this  Bessel 
function,  a  numerical  study  was  made  of  t'  .  From  the  power 
series  definition  of  Kq  ,  and  a  term-by-term  integration  of 
Equation  (3-29)  ve  obtain 
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where  y  =  .5772156649...,  (Euler's  constant) 


By  use  of  recursion  relations  between  the  coefficients,  the 
series  is  readily  programmable,  with  high  accuracy,  even  on  a 
pocket  calculator,  for  values  of  z  as  large  as  12.  However, 
for  z  greater  than  about  8,  higher  accuracy  may  be  achieved 
by  an  asymptotic  series.  By  expanding  the  exponential  in  a 
power  series  and  integrating,  the  following  expressions  are 
obtained 

00 

P  -  /X  exp(-ab£  (-l)j  dj  a‘i(i+%)  ,  (3-32) 
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This  series  is  "semi -convergent"  and  must  be  truncated  when  the 
magnitude  of  the  terms  no  longer  decreases.  For  a^  greater 
than  about  8,  it  is  capable  of  providing  3-plac.e  accuracy  for  F' 


deteriorates;  Lor  a  ft;  1  ,  i  L  is  eiiecLiveiy  worthless. 
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The  f tine t ion  t'  =  1-  ~  J  KQ(x)dx  is  plotted  in  Figure  3- 

o  _ 

in  Cartesian  coordinates  and  the  curve-of-growth  -*.n  t'  is  plotte 
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Figure  3-2.  Curve- of -Growth  -Jln(t')  versus  a  for 

f.  ,  1  -  i  /*\.0(x)dx  . 

o 


3.2.2  Random  band  models:  The  mean  trasnmissivity  for  a 
given  band  model  may  be  written  alternatively 
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(3-34) 


in  which  f(k)  is  the  probability  density  function  of  the 
absorption  coefficient  k  .  Because  the  properties  of  random 
band  models  without  scattering  are  relatively  simply  derivable 
from  those  of  single  lines,  explicit  expressions  for  such 
probability  densities  have  not  been  required  for  study  of  the 
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gaseouB  band  models.  However,  because  of  the  exponential  form 
of  the  monochromatic  transmissivity,  the  probability  density 
function  in  the  expression 

oo 

t(u)  "  y*  f(k)e~ku  dk  (3-35) 

o 

can  be  seen  tc  be  the  inverse  Laplace  transform  of  the  mean 
transmissivity  function.  A  study  of  the  inverse  Laplace  trans¬ 
forms  of  known  transmissivity  functions  has  been  made  by 
Domoto  (1974),  who  has  found  that  for  the  exponential-tailed 
S'1  model  a  relatively  simple  functional  form  exists,  namely, 
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corresponding  to  the  transmissivity  function 
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We  consider  here  the  case  of  no  particulate  absorption,  for  which 
the  effects  of  the  gaseous  spectral  structure  will  be  the  greatest. 
The  function  t’  previously  defined  as 
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ku)  can  be  rewritten 
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Direct  integration  is  rendered  difficult  by  the  presence  of 
the  radical  in  the  exponential  term.  Before  studying  the  general 
behavior  of  t ’ ,  some  limiting  cases  (small  t8.  large  rfl)  are  examined. 

a)  Small  t_ 

Q 

For  small  values  of  xg  (rg  <<  xa)  ,  expansion  of  the  radical 
yields 
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in  which  terms  in  etc.  have  been  suppressed.  Substitution 

produces 
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Since  the  integral  is  now  in  the  standard  form  with  new  constants, 
we  have 
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or 
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To  the  same  degree  of  approximation  we  may  write 
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b )  Large  t 


In  the  limit  of  large  t  the  magnitude  of  the  argument 

s 

of  the  exponential  term  in  the  integral  becomes  large,  and  the 
minimum  progressively  less  shallow,  and  occurring  at  smaller 


values  of  x 
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the  argument  of  the  exponential  terra  becomes 
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The  magnitude  of  this  argument  has  a  minimum  at  approximately 
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so  that  for  large  x  ,  t  min  <<  t„ 
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In  the  region  of  this  minimum,  the  argument  may  be  expanded: 
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the  expression  for  t'  then  becomes 
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By  recourse  to  the  mean  value  theorem,  noting  the  narrowness 
of  the  maximum  of  the  exponential  in  the  neighborhood  of 

ttW  _  ^0/9 

t  ,  the  multiplicative  factor  of  t  “  '  may  be  taken 
outside  the  integral  and  evaluated  at 
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The  remaining  integration  may  be  performed  by  letting  the  lower 
bound  -+-<»: 


t'  exp  -3(%axaig) 


(3-54) 


This  functional  dependence  provides  an  interesting  analogy  to 
the  square-root  region  for  pure  gaseous  absorption  in  which 


t  exp  I  -2 .(a ia)‘ 
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a)  Arbitrary  x^ 

& 

Numerical  studies  have  been  carried  out  on  the  integral 
representation  of  the  function  IT',  given  by  Eq.  (3-39).  Some 
results  of  these  calculations  are  shown  in  Figures  3-3  to  3-5, 
in  which  t'  is  plotted  against  x  ,  the  mean  absorption 

a 

optical  depth,  for  various  values  of  x  ,  the  scattering 

optical  depth  and  a  ,  the  fine  structure  coefficient.  The 

linear  scale  for  x'  does  not  show  the  precise  asymptotic  behavior 

at  x*  very  close  to  0  or  1,  but  does  indicate  a  smooth,  well- 

behaved  dependence  on  x  ,  encouraging  confidence  in  the  approach 

s 

of  seeking  an  analytical  or  semi -empirical  functional  represen¬ 
tation. 


i-t'  r 


Figure  3-3.  Plot  of  t'  vs.  x*  for  a  =  tt/1000  and 
x  ranging  from  0  to  10^. 

S 
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Figure  3-5.  Plot  of  P  vs.  ?a  for  a  =  tt/10  and 
xs  ranging  from  0  to  102. 


O  rrj 


The  asymptotic  region  of  small  t'  was  examined  in  detail. 


The  conclusion  was  reached  that  the  particular  asymptotic  form 
developed  above  | see  Eq.  (3-54)]  is  of  major  practical  significance 
A  plot  of  the  ratio  t:' /exp  [-3  (%a7  t  J*5 ]  is  presented  in  Figure 
3-6.  This  ratio  can  be  seen  to  approach  the  predicted  theoretical 

JL 

value  of  (4/3)^  =  1.1547  closely  as  a-*0  even  for  small  values  of 


Figure  3-6. 


Plot  of  t  /exp[-3(%aT  x  )  ]  vs .  at  t  showing 

asymptotic  behavior  for  "large"  t  (ranging  from 

3  S  — 

1  to  10  ) ,  for  various  values  of  a  and  ax  . 

_  a 

The  range  of  values  of  ax  x  corresponds  to  a 

as 

range  of  t'  from  about  0.7  at  the  left  to 
essentially  zero  at  the  right. 
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af  t  ,  when  t  is  sufficiently  large.  Empirically,  for  a 
as  s 

<0.05  and  t X  10ax  ,  the  asymptotic  approach  is  quite  close 

^  2 

(within  about  10%,  even  for  a?aTg  as  small  as  10  ,  corres¬ 

ponding  to  transmissivities  as  large  as  0.6).  For  larger  values 
of  a,  an  asymptotic  value  is  approached,  whose  magnitude 
depends  on  the  value  of  a.  The  particular  functional  dependence 
could  be  determined  theoretically  or  empirically  (e.g.,  an 
asymptotic  value  of  (4/3)^exp (2a)  provides  an  excellent  fit  over 
the  range  of  values  of  a  considered  here) . 


3.2.3  Empirical  formulations  for  scattering  band  models. 

A  study  was  made  of  possible  empirical  forms  in  the  scattering 
band  model  representation  based  on  the  previous  theoretical 
studies  and  numerical  calculations.  No  exact  explicit  expres¬ 
sions  were  derived  for  the  general  case,  but  various  limiting 
cases  were  studied,  and  any  empirical  form  considered  should 
have  the  appropriate  theoretical  asymptotic  behavior.  The 
numerical  calculations  then  serve  to  select  the  best  representa¬ 
tion  for  intermediate  values  of  the  parameters. 

Specifically,  the  function  previously  defined  at  t' 
(analogous  to  the  mean  transmissivity  function  t  for  the 
no-scattering  limit)  is  desired  to  have  the  following  limits: 
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A)  No  scattering  (t  ”0)  •  t'  ■*  t  =>  exp  |-f(t  ,  a)|  where  f 
is  the  appropriate  gaseous  band  model  curve  of  growth, 
in  which  a  is  the  fine- structure  parameters  and  ?a 
the  mean  absorption  optical  depth  (»  ku  where  k  is  the 
gaseous  absorption  coefficient  and  ux  the  amount  of  the 
absorbing  gas,  in  reciprocal  units).  The  function  f 
may  include  particulate  absorption,  i.e.,  be  replaced  by 
f(t  ,  a)  +  t  where  f  is  the  pure  gas  curve  of  growth  and 
tq  is  a  spectrally  constant  absorption  optical  depth 
(gaseous  and/or  particulate). 

B)  Weak  scattering  (t  4:4  7  ) 

t'  ^  exp  -f(?a,  a)  - t g 

in  which  f  may  include  particulate  absorption  as  in  (A) . 

C)  Constant  absorption  coefficient  (gaseous  and/or  particulate). 

(a  -*•  00  ;  f(7r  ,  a)  -*•  x  for  any  band  model) 

3  3 

t'  exp  -(t  z  +  2t  t  )^  (this  is  a  func- 
c  a  as 

tional  requirement  of  either  the  two- flux  or  six- flux 
treatment  for  constant  absorption  coefficient). 

D)  Strong- line,  large  scattering  limit  (no  particulate  absorp¬ 
tion) 

(?a/a  large,  xQ  large) 


for  a  random  band  model,  with  cp'-l  and  c?''-2.4.  This 

important  limiting  behavior  was  derived  in  the  previous 

study  and  is  analogous  to  the  strong-line  limit  for 

\ 

gaseous  absorption  t  ►  exp  -2(i  a)  J 

a 

An  empirical  form  satisfying  (A) ,  (B) ,  and  (C)  is 


t' 


exp 


a) 


(3-55) 


This  is  the  form  utilized  intuitively  in  the  Vanderbilt-Slack 
program.  However,  in  limiting  case  (D)  it  approaches 
t'  +  exp (-27a%  a^  ts*)  . 


Alternatively,  a  form  satisfying  (A),  (B) ,  and  (D)  is 


t'  =  exp  -f(ia,a)  jl+ 


f(Ta-a)l 


jl/3 


(3-56) 


Although  this  form  provides  a  better  fit  to  the  numerical  values 
calculated  for  intermediate  ranges  of  the  parameters,  it  does  not 
have  the  proper  form  to  match  limit  (C) ;  instead 

—  —  —  I/O 

t'  -►  exp  -(x  3+  3xa2is)  •  For  this  reason,  a  more  elaborate 
form  containing  all  four  limiting  cases  above  is  sought. 


A  convenient  form  having  these  required  limits  is 


.gLui  .jin i  . .  . . 


where  b,  c  and  n  are  empirical  constants.  Numerical  evaluation 
of  this  equation  has  shown  significantly  improved  agreement 
relative  to  the  previous  forms.  The  curve  of  growth  f(?  ,a) 
may  include  a  particulate  absorption  optical  depth. 

Figure  3-7,  -9,  and  -11  show  comparisons  of  the  exact  values 

of  t'  previously  calculated  for  ir/103  <  a  <  ir/10  (solid  lines) 

with  the  approximation  given  by  Equation  (3-55)  (dashed  line) . 

This  approximation  is  seen  to  be  especially  poor  for  large 

t  and  small  a  (i.e.,  when  condition  (D)  should  be  met), 
s 

Figure  3-8,  -10  and  -12  show  comparisons  of  the  same  values 
of  t'  with  the  empirical  form  given  by  Equation  (3-57)  for  the 
case  c“l,  n=l,  b=l/8.  The  improved  match  overall  is  evident. 

No  effort  was  made  at  optimization  of  the  parameters  b,  c,  and 
n.  It  is  anticiapted  that  an  optimized  version  of  this  form 
(or  a  very  similar  one)  will  be  used  in  the  basic  formulation 
of  the  SIRRM  band  model  representation. 

All  the  results  shown  here  have  been  obtained  for  the 
exponential-tailed  1/S  random  band  model.  The  expression  for 
the  absorption  coefficient  probability  distribution  function 
for  the  Goody  model  is  so  mathematically  intractable  as  to  make 
it  extremely  unattractive  for  similar  theoretical  or  numerical 
studies.  However,  the  conclusions  reached  here  regarding  proper 
asymptotic  and  approximate  forms  are  certainly  also  applicable 
to  the  Goody  model  (whose  curve  of  growth  for  the  non-scattering 
case  is  only  slightly  different  from  that  of  the  exponential¬ 
tailed  1/S  model). 
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-  n/iOOO 


!  10”1  1  10 

10J  103 

T. 

Figure  3-7.  Plot  of  t'  vs.  ia 

fo: 

r  a  =  tt/1000  and 

t  g  ranging  from  0 

to 

10 4  compared  with 

approximation  (1) 

Figure  3-8 


Plot  of  t'  vs.  t  for  a  -  tt/1000  and 

a 

t  ranging  from  0  to  10“  compared  with 

S 

approximation  (3). 
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vs.  tg  for  a  =  tt/100 
from  0  to  10 3  compar 
on  (3). 


t  -iotn 


5  x. 


a  “  «/10 


Figure  3-11.  Plot  of  t  vs.  t  for  a  =  ir/10  and 

a 

tg  ranging  from  0  to  10 2  compared  with 
approximation  (1). 


Figure  3-12. 


Plot  of  t'  vs.  t  for  a  *  tt/10  and 

t  ranging  from  0  to  10 2  competed  with 
s 

approximation  (3), 
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The  preceding  development  is  for  single  band  models  (with 
or  without  particle  or  other  gaseous  absorption).  In  practice, 
it  is  anticipated  that  methods  of  combining  band  model  parameters 
will  be  utilized  which  are  analogous  to  those  in  current  use 
in  gaseous  radiative  transfer  codes  (e.g.,  the  Curtis-Godson 
approximation) . 
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OPTICAL  PARAMETERS  OF  PARTICLES 
4.1.  Introduction 

The  treatment  of  the  gas/particles  radiative  transfer  requires 
as  inputs  the  volume. extinction  coefficients  and  angular  scat¬ 
tering  phase  functions.  These  quantities  can  be  derived  from  the 
primary  optical  parameters  by  use  of  the  Mie  theory,  which  des¬ 
cribes  the  scattering  and  absorption  properties  of  a  homogeneous 
sphere  whose  optical  properties  are  given  by  the  complex  index 
of  refraction.  Although  the  Mie  theory  is  based  upon  the  ideal¬ 
ization  of  homogeneous  spherical  particles,  it  can  be  considered 
as  a  very  reasonable  first  approximation  since  real  plume  particles 
were  often  found  to  be  solid  spheres  (Rochelle,  1967).  In 
addition,  higher  order  account  of  particle  non-uniformity  and 
non-sphericity  is  prohibitively  complex  (Kerker,  1969)  and  not 
warranted  by  the  present  state  of  knowledge.  Exact  solutions 
of  the  wave  equations  are  obtained  by  the  matching  of  boundary 
conditions  at  the  surface  of  the  sphere.  The  solutions  are  ex¬ 
pressed  as  infinite  series  involving  Legendre  polynomials  and 
Picatti-Bessel  functions.  Evaluation  of  the  formal  series 
solution  proves  computationally  difficult  for  some  ranges  of 
the  parameters  and  considerable  effort  has  been  expended  in  the 
past  on  the  development  of  efficient  and  accurate  computational 
methods  of  evaluating  these  solutions.  The  most  recent  treat¬ 
ment  is  by  Dave  (1972)  and  it  is  his  approach  that  is  being  used 
in  the  present  program.  A  brief  discussion  together  with  the  ex¬ 
tension  from  single  particles  to  particle  distributions  is  given 
in  Section  4.2. 
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A  preliminary  listing  of  the  particulate  species  of  interest 
is  given  in  Table  4-1  together  with  an  indication  of  the  avails- 

t 

bility  of  the  index  of  refraction  (ni )  and  the  absorption 
index  (n2)  measured  in  a  given  temperature  interval.  Only  four 
out  of  nine  particulate  species  are  sufficiently  known,  and 
these  will  have  to  serve  as  representatives  until  more  informa¬ 
tion  about  the  remaining  species  becomes  available. 


TABLE  4-1.  LIST  OF  EXHAUST  PARTICULATE  SPECIES 
AND  THE  AVAILABILITY  OF  ^  and  n2 . 


PARTICLE 

n1(T,A) 

n2 (T, A) 

SOURCE 

A£203 

(s) 

300-1800K 

300-1800 

f 

Aerospace 

t 

(s) 

300K 

300-3000 

Aerospace/PRA  Estimate 

MgO 

(s) 

300-2225K 

300-2225 

Aerospace 

Zr02 

(s) 

300-600K 

300-600 

Aerospace 

ZrC 

(s) 

BN 

(s) 

300K(Na-D) 

Kordes  (1960) 

B2°3 

00 

300K(Na-D) 

Handbook  Chem. Phys . for  (s) 

BeO 

(s) 

300K(Na-D) 

Handbook  Chem . Phys . 

Be3N2 

Representative  results  of  volume  coefficients  and  angular 
scattering  phase  functions  for  and  carbon  particle  distri¬ 

butions  are  given  in  Section  4.3.  The  results  indicate  that  the 
large  variations  in  coefficients  and  phase  functions  typically 
observed  as  a  function  of  wavelength  and  scattering  angle  for 
single  particles  are  greatly  reduced  when  realistic  size  dis¬ 
tributions  are  introduced.  In  addition,  the  influence  of  the 
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A  preliminary  listing  of  the  particulate  species  of  interest 
is  given  in  Table  4-1  together  with  an  indication  of  the  availa- 

i 

bility  of  the  index  of  refraction  (n! )  and  the  absorption 
index  (n2)  measured  in  a  given  temperature  interval.  Only  four 
out  of  nine  particulate  species  are  sufficiently  known,  and 
these  will  have  to  serve  as  representatives  until  more  informa¬ 
tion  about  the  remaining  species  becomes  available. 


TABLE  4-1.  LIST  OF  EXHAUST  PARTICULATE  SPECIES 
AND  THE  AVAILABILITY  OF  ^  and  n2 . 


PARTICLE 

nL (T, A) 

n2(T,A) 

SOURCE 

A£^0q 

(s) 

300-1800K 

300-1800 

Aerospace 

£ 

(s) 

300K 

300-3000 

Aerospace/PRA  Estimate 

MgO 

(s) 

300-2225K 

300-2225 

Aerospace 

Zr02 

(s) 

300-600K 

300-600 

Aerospace 

ZrC 

(8) 

BN 

(8) 

300K(Na-D) 

Kordes  (1960) 

U) 

300K(Na-D) 

Handbook  Chem. Phys . for  (s) 

BeO 

(s) 

300K(Na-D) 

Handbook  Chem. Phys. 

Be3N2 

Representative  results  of  volume  coefficients  and  angular 
scattering  phase  functions  for  Ai^O^  and  carbon  particle  distri¬ 
butions  are  given  in  Section  4.3.  The  results  indicate  that  the 
large  variations  in  coefficients  and  phase  functions  typically 
observed  as  a  function  of  wavelength  and  scattering  angle  for 
single  particles  are  greatly  reduced  when  realistic  size  dis¬ 
tributions  are  introduced.  In  addition,  the  influence  of  the 
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4 

OPTICAL  PARAMETERS  OF  PARTICLES 
4.1.  Introduction 

The  treatment  of  the  gas /particles  radiative  transfer  requires 
as  inputs  the  volume  extinction  coefficients  and  angular  scat¬ 
tering  phase  functions.  These  quantities  can  be  derived  from  the 
primary  optical  parameters  by  use  of  the  Mie  theory,  which  des¬ 
cribes  the  scattering  and  absorption  properties  of  a  homogeneous 
sphere  whose  optical  properties  are  given  by  the  complex  index 
of  refraction.  Although  the  Mie  theory  is  based  upon  the  ideal¬ 
ization  of  homogeneous  spherical  particles,  it  can  be  considered 
as  a  very  reasonable  first  approximation  since  real  plume  particles 
were  often  found  to  be  solid  spheres  (Rochelle,  1967).  In 
addition,  higher  order  account  of  particle  non-uniformity  and 
non-sphericity  is  prohibitively  complex  (Kerker,  1969)  and  not 
warranted  by  the  present  state  of  knowledge.  Exact  solutions 
of  the  wave  equations  are  obtained  by  the  matching  of  boundary 
conditions  at  the  surface  of  the  sphere.  The  solutions  are  ex¬ 
pressed  as  infinite  series  involving  Legendre  polynomials  and 
Ricatti-Bessel  functions.  Evaluation  of  the  formal  series 
solution  proves  computationally  difficult  for  some  ranges  of 
the  parameters  and  considerable  effort  has  been  expended  in  the 
past  on  the  development  of  efficient  and  accurate  computational 
methods  of  evaluating  these  solutions.  The  most  recent  treat¬ 
ment  is  by  Dave  (1972)  and  it  is  his  approach  that  is  being  used 
in  the  present  program.  A  brief  discussion  together  with  the  ex¬ 
tension  from  single  particles  to  particle  distributions  is  given 
in  Section  4.2. 
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temperature  on  the  coefficients  and  phase  functions  for  solid 
carbon  is  minimal . 

In  the  last  section,  we  recommend  a  scheme  for  the  Data 
File  that  combines  flexibility  and  minimum  requirements  in  the 
operation  of  the  SIRRM  code. 

4 . 2  Formalism  of  the  Hie  Scattering 

4.2.1  Single  particle.  In  the  Mie  theory,  the  primary 
quantities : 

Si (a,n,0) 

S2 (a,n,0) 

<W(a’n> 

%ca(a'n> 

are  calculated  (Van  de  Hulst,1957)  where  S  and  S  are  amplitude 

1  2 

functions  related  to  the  planes  of  polarization,  Q  and  Q 

cXt  SCa 

are  the  extinction  and  scattering  efficiencies,  a  is  the  size 
parameter  (a=2Trr/A,  with  radius  r  of  scattering  sphere  and  wave¬ 
length  A),  n  is  the  refractive  index  (n=n(A)  =  nj  (A) -in2 (A) )  and 
0  is  the  angle  between  the  directions  of  incoming  and  scattered 
radiation.  After  the  primary  quantities  are  calculated,  the 
intensities  of  the  scattered  radiation  and  their  degrees  of 
polarization  can  be  determined.  The  calculation  is  simplified 
by  repeated  use  of  recursion  relationships  and  efficient 
computing  methods  have  been  developed  (Deirmendj ian,  1969;  Dave 
1969,  1970A,  1970B) .  In  the  following,  the  approach  taken  by 
Dave  (1972)  is  briefly  described. 
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The  intensity  scattered  by  a  single  spherical  particle 
per  unit  solid  angle  for  unity  incident  intensity  is  given  by 

I(a,n, 0)  ("r2)  Qs(a,n)  M(a,n,0)  (4-1) 

where  M(a,n,e)  is  the  normalized  scattering  phase  function,  which 
is  given  by 

n 

M(a,n, e)  Afc(a,n)  Pk_1(cos  0)  (4-2) 


where  Pkl(cos  0)  are  the  Legendre  polynomials.  The  coefficients 
of  the  Legendre  series  are  given  by 


Ak(a,n) 


2 

a'2Qs  (a,n) 


r  a) 

Lk  (a ,n) 


(2) 

+  Lk  (a  ,n) 

. 


where  the  functions  1^(1) . (2)  (a,n)  are  computed  from  the  values 
of  the  complex  Hie  amplitude  an(a,n)  and  bn(a,n).  The  extinction 
and  scattering  efficiencies  are  given  by 

W-n>  -  5  £1<2n+1>  Re{*n  +  V 

Wa’n>  '  j;  £1<Zn+1)  +  'bnl!) 

The  amplitudes  an  and  bn  are  expressible  in  terms  of  the  spherical 
Bessel  functions  of  the  first  and  second  kind.  Using  the  method 
of  logarithmic  derivative  function  An(a,n)  of  the  Bessel  functions, 
the  numerical  computations  are  facilitated.  However,  difficulties 
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exist  in  computing  A^Ca.n)  under  certain  conditions  (Deirmendj ian , 
1969).  Kattawar  and  Plass  (1967)  found  that  these  difficulties 
can  be  avoided  by  using  a  downward  recurrence  scheme  for  A^Ca.n) 
which  is  basically  stable.  Dave  (1972)  has  adopted  this  scheme 
in  his  approach. 

4.2.2  Particle  size  distribution.  Since  rocket  exhaust 
plumes  do  not,  in  general,  consist  of  particles  with  only  one 
fixed  size,  it  is  necessary  to  investigate  the  characteristics 
of  size  distributions.  These  characteristics  are  described  by 
the  volume  extinction  coefficient  6  ,  the  volume  scattering 

coefficient  8sca>  the  volume  absorption  coefficient  3a^s  and 
the  normalized  scattering  phase  function  M(0) .  These  quantities 
are  dependent  on  n(A)  and  various  parameters  describing  dif¬ 
ferent  size  distributions.  They  can  be  expressed  as 

K 

M(0)  -  ]£  Ak(r-)  Pk-l^cos  °')  (4-4) 

k=l  J 

and 

rj 

32(rj)  =  it  f  Qz(a,n)r2  N(r)dr  (4-5) 

rmin 
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where  z  stands  for  ext,  sea  or  abs.  The  normalized  Legendre 
coefficients  become 


A  (r  )  =  - 

K  J  7T  8 


scav  j 


j 

-  f  Li(a*n)  N(r)dr 

(r.)  J  k 


(4-6) 


The  phase  function  is  normalized  according  to 


yH(e)  sin0d9d({>  =  4n 


The  distribution  N(r)  is  given  by 


N(r)  =  Nf (r) 


(4-7) 


where  N  is  total  number  of  particles  per  unit  volume  and 
00 

I"  f(r)dr  =  1.  Dave  (1972)  has  made  provisions  to  use  any 
) 

one  of  three  different  size  distribution  functions.  We  are 
interested  in  the  discontinuous  power  law  type 


f (r)  -  C 


for  r  .  <  r  <  r 

ml  ^  ^  « 


f(r>  =  C(rnl/r)''+1  for  rm  <  r  «  rmax 


(4-8) 
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and  the  modified  gamma  function 

f(r)  »  a  rVbrYfor  rmln  <  r  «  rBax  (4-9) 

It  is  found  that  the  modified  gamma  function  is  representative 
of  9ize  distributions  encountered  in  exhaust  plumes  with 
A£20j  particles  (Dawbarn,  1978).  However,  many  different  size 
distributions  have  been  encountered  which  sometimes  depend  upon 
the  method  of  collection  (see  Figure  4-1) .  The  coefficients 
a,  a,  b  and  y  can  be  adjusted  to  represent  any  of  the  measured 
or  calculated  size  distributions. 


Particle  Sire,  pm 


Figure  4-1.  Two  A£203  particle  size  distributions  of  Titan  IIIc 

using  two  different  collection  methods  (Dawbarn,  1978) 
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The  parameters  for  the  modified  gamma  function  are  not  in¬ 
dependent  of  each  other  and  they  have  to  be  adjusted  in  a  certain 
manner  to  obtain  a  given  size  distribution  that  has  its  maximum  at 


radius  rm  .  This  interdependence  was  investigated.  The  four 
parameters  may  be  reduced  to  two  by  utilizing  the  conditions 


=  0  for  r  =  r 


(4-10) 


and 

oo 

/  f(r)  =  1 

c 

These  conditions  lead  to 

a/b  =  Y  rmY  (4-11) 

1  *  a  y  ^  b  ^r(()  (4-12) 

where  £  =  —  and  r(0  is  the  gamma  function.  At  r  =  rm> 

f(r)  must  assume  a  value  between  zero  and  unity  (as  determined  by 
experiment).  Introducing  Equations  (4-11)  and  (4-12)  into  (4-10) 
and  solving  for  r  results  in 


(4-13) 


Assuming  upper  and  lower  values  for  f(rm)  .  Equation  (4- 13)  has 
been  solved  for  f (r^l  and  f(rm)^.2  using  different  values 
of  a  and  y  •  The  solutions  are  shown  graphically  in  Figure  4-2. 
A  family  of  curves  is  made  up  for  a  given  a  ,  the  upper  and 


4-0 


lower  limits  of  which  are  given  by  f(rm)'vl  and  f(rm)'v.2  . 

Thus,  pairs  of  a  and  y  can  be  chosen  to  describe  any  size 

distribution  for  r  between  0.01  and  10  |.tm .  The  choice  is 

m 

dictated  by  the  known  (measured)  standard  deviation  of  the  dis 
tribution  function.  The  higher  values  of  a  and  y  produce 
smaller  standard  deviations . 


Figure  4-2.  Families  of  Curves  for  which  0.2*f(rm)'  l 
at  a  given  r  as  determined  by  pairs 
of  y  and  a. 
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Following  a  similar  argument  as  was  presented  above  for  the 
modified  gamma  function,  the  adjustable  parameters  in  the  dis¬ 
continuous  power  law  (Eq.  4-8)  are  determined  through 


f(r) 


dr  +  C  r_ 


ou 

/  r 


(v+1) 


1/C-r 


Based  on  a  review  of  existing  experimental  data  (Rochelle,  1967), 

we  have  chosen  r  =>0.03  pm  for  solid  carbon  particles  and  C  =  1. 

m 

With  these  assumptions,  then  v^O.  The  size  distribution  is 
shown  in  Figure  4-3.  It  resembles  that  used  by  Stull  and  Plass 
(1960)  which  was  based  on  data  by  Tesner  (1959).  Included  in  the 
same  figure  is  also  the  size  distribution  based  on  Equation  (4-9) 


■■••fttfaaal 


Figure  4-3. 


Size  Distribution  of  Solid  Carbon  Particles  in  Rocket 
Exhaust,  based  on  the  Two  Distribution  Functions  Des¬ 
cribed  in  the  Text.  Distribution  used  by  Stull  and 
Plass  (1960)  is  indicated  by  (+) . 


4 . 3  Optical  Paramo  tars 


A  set  of  representative  values  of  ni  and  n2  were  collected 
by  the  Aerospace  Corporation  for  AH203,  C,  MgO,  and  Zr02 
(Whitson,  1975).  Additional  data  for  these  particles  were  pre¬ 
sented  by  Dowling  and  Randall  (1977)  in  part  based  on  experiments 
conducted  during  the  "Particle  Optical  Properties  Measurement 
(POPM)  Program".  Our  effort  was  concentrated  in  comparing  these 
data  and  study  the  influence  of  higher  temperatures. 

Carbon 

As  an  example  of  the  wide  differences  of  measured  refractive 
index  nt  ,  and  absorption  n2 ,  Dowling  and  Randall  have  compared 
data  on  polycrystalline  and  pyrolytic  graphite,  pressed  carbon  soot, 
and  a  microcrystalline  polished  slab  of  graphite.  The  results 
for  n2  are  shown  in  Figure  4-4.  We  have  added  the  NASA  Handbook 


Figure  4-4.  The  Absorption  Index  of  Carbon  Materials  in 
the  Wavenumber  Region  400  cm-1  to  2200  cm-!. 
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data  for  room  temperature.  The  "representative"  values  chosen 
by  Whitson  are  those  for  the  polycrystalline  graphite. 

In  addition  to  room  temperature  data,  the  NASA  Handbook 
lists  high  temperature  data  for  the  absorption  coefficient  of  small 
carbon  particles.  These  data  were  based  in  part  on  theoretical 
work  and  in  part  on  measurements.  We  have  extrapolated  those 
data  to  3000  K  and  to  A  «  25  pm,  using  the  relationship 


to  convert  the  absorption  coefficients  kc  (cm?/gm)  into  absorption 
indices,  where  pQ  =  2g/cm3  was  chosen  as  a  representative 
value  for  the  density  of  carbon.  The  results  are  plotted  in 
Figure  4-5.  We  recommend  to  use  these  values  for  n 2  until  more 
experimental  data  become  available.  For  large  particles  compared 
with  A,  the  computed  particle  emissivity  approaches  the  emissivity 
of  a  plane  bulk  sample.  Dowling  and  Randall  have  calculated 
en  -  .71  for  A  <  25  pm,  choosing  nt  =  3  and  n2  =  1 .  Using  the 
same  formula,  we  have  calculated  en  parametrically  for  0  «  n^  4 
and  0.5  <  n2  <  1.5.  The  results  are  shown  in  Figure  4-6. 

Measured  spectral  smissivities  as  reported  by  Whitson  range  from 
.8  to  1.0  for  different  samples  of  graphite  (except  for  pyrolytic, 
c-face)  in  the  temperature  regime  from  1100  to  2300  K,  pointing 
to  a  somewhat  lower  value  of  n i  at  higher  temperatures. 


The  only  high  temperature  data  for  n  listed  by  Whitson 
are  due  to  N.  K.  Krascella  (1966),  and  are  plotted  in  Figure  4-7 
together  with  Whitson's  "representative"  values  at  room  tempera¬ 
ture  and  our  extrapolation  ’to  longer  wavelengths. 


Figure  4-7.  nx (Carbon)  vs.  A  for  T  »  300  and  2250K, 
based  on  the  Compilation  by  Whitson 
and  our  Extrapolation. 

A£,03 

Experimental  reflectivity  measurements  in  the  POPM  program 
at  300  and  678  K  were  used  to  determine  n1  and  n2  of  A£703 
(Dowling  and  Randall,  1977).  The  results  are  shown  in  Figure 
4-8  and  4-9  and  they  indicate  that  nx  and  n2  are  reduced  at 
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300  K 


Figure  4-8.  The  Refractive  Index  (n,)  of  A£203  in  the 

Region  400  cm"1-  to  2200  cm-^  as  Determined  from 
a  Kramers -Kronig  Analysis  of  Reflectivity  Data 
taken  at  T=300  and  678  K  (Dowling  and  Randall,  1977). 


Figure  4-9.  The  Absorption  Index  (n2)  of  A£,,03  in  the 
400  cm"!  to  2200  cm"1  as  Determined  from  a 
Kramers -Kronig  Analysis  of  Reflectivity  Data 
taken  at  T-300  and  678  K  (Dowling  and  Randall,  1977) 
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higher  temperatures  in  the  resonance  region.  A  comparison  with 
the  "representative"  values  by  Whitson  shows  good  agreement  for  the 
room  temperature  data.  Also  the  qualitative  behavior  at  higher 
temperatures  in  the  resonance  region  has  been  observed  by  B.  Piriou 
(1966).  The  data  as  presented  by  Whitson  are  shown  in  Figures  4-10 
and  4-11. 


Figure  4-10.  nj(A£203)  as  a  Function  of  X 

for  Two  Temperattires  as  Pre¬ 
sented  by  Whitson  (1975) . 
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Figure  4-11.  n2(A£203)  as  a  Function  of  A 

for  Two  Temperatures  as  Pre¬ 
sented  by  Whitson  (1975) . 


MgO 


M.  E.  Whitson  (1975)  has  listed  experimental  data  of 
ni  and  nz  taken  by  several  workers.  All  of  these  are  in  general 
agreement  and  the  "representative1'  values  at  room  temperature 
are  based  on  these  values.  There  appears  to  be  little  or  no 
dependency  on  temperature,  as  the  data  by  B.  Piriou  (1966) 
indicate  (see  Figure  4-12  and  4-13).  The  experimental  study 
and  subsequent  data  analysis  using  classical  dispersion  theory 
by  Dowling  and  Randall  (1977)  confirm  these  data.  A  comparison 
of  ni  and  n2  with  Whitson's  representative  room  temperature  data 
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Figure  4-12.  n  of  MgO  vs.  X  for  Three  Different  Temperatures 


show  good  agreement  (see  Figures  4-14  and  4-15).  In  addition, 
Dowling  and  Randall's  data  at  elevated  temperature  (573  K)  show 
little  or  no  temperature  dependency  (see  Figures  4-16  and  4-17) . 

Whitson's  report  does  not  contain  any  data  of  m  and  only 
a  limited  set  of  n2  at  room  temperature  between  5  and  9  um. 

Thus,  the  data  obtained  by  Dowling  and  Randall  filled  an  important 
gap.  The  reflectivity  between  400  and  2200  cm”'  was  measured 
at  room  temperature  and  573°K.  The  spectral  structure  is  similar 
to  that  reported  by  Whitson,  albeit  lower  by  about  207..  The 
results  for  ni  and  112  are  shown  in  Figures  4-18  and  4-19, 
respectively . 

4 . 4  Numerical  Results  for  AftpOo  and  Carbon 

The  volume  extinction,  scattering  and  absorption  coefficients 
and  the  scattering  phase  function  have  been  calculated  for  AA203 
and  carbon  with  a  code  developed  by  Dave  (1972).  These  two 
species  are  the  most  important  particulate  constituents  in 
present  missile  exhausts  and  possess  very  different  optical  and 
geometrical  properties.  Thus,  the  results  give  insight  into 
the  functional  dependencies  of  8  and  M  on  size  distributions, 
wavelength  and  temperature.  For  A&2O3,  the  modified  gamma 
function  was  used  for  representing  the  size  distribution  2  as 
given  in  Figure  4-1,  viz. 

N(r)  -  N  .0551  y  e" * 004r4 
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(n2)  of  MgO  from  400  cm 
red  by  Classical  Dispersion 


~  3 

The  total  number  of  particles  chosen  was  10  cm  which  cor¬ 
responds  to  a  mole  fraction  of  about  307.  when  the  mean  radius 
is  5  pm.  The  min  and  max  radii  chosen  were  1  and  9  ym, 
respectively.  The  input  for  affi^n>  amax,  nj  and  n2  as  a  function 
of  X  for  a  temperature  of  1773°K  are  given  in  Table  4-2. 

TABLE  4-2.  INPUT  PARAMETERS  USED  FOR  A*203 


X(pm) 

a  min 

a 

max 

ni 

n2 

1 

6.29 

56.55 

2 

1.5-3 

3 

2.09 

18.85 

2 

2. 0‘3 

5 

1.26 

11.31 

2 

6 . 0-3 

7 

.90 

8.08 

1.5 

1.8"2 

9 

.70 

6.28 

1.2 

8. 0'2 

11 

.57 

5.14 

.5 

.3 

13 

.48 

4.35 

.3 

1.0 

15 

.42 

3.77 

.3 

2.5 

17 

.37 

3.  33 

.5 

6.0 

19 

.33 

2.93 

4.0 

.6 

21 

.30 

2.69 

3.5 

8 

23 

.27 

2.46 

3.0 

6.0 

25 

.25 

2.26 

5.0 

1.0 

The  results  for  the  volume  extinction,  scattering  and  ab¬ 
sorption  coefficients  and  albedo  (ui  =  6„„a/BovJ  as  a  function 
of  X  (in  steps  of  2  ym)  are  listed  in  Table  4-3  and  the  angular 
scattering  functions  M(9)  as  functions  of  0  are  shown  In  Figures 
4-20  through  4-32.  It  can  be  seen  that  the  scattering  becomes 
more  and  more  isotropic  as  the  wavelength  increases.  This  is  also 
shown  in  Figure  4-33,  where  the  ratio  of  the  side  to  forward 
scattering  component  is  plotted  vs.  A. 


TABLE  4-3.  RmuUs  for  8„t .  B,c,.  8abs 

and  u  as  a  Function  of  A  for 
Afc2Os  Particles  in  Rocket 
Plume  Exhaust  at  1773K. 


1<V«> 

*•» 

®*b» 

fct 

1 

2.182 

1.961 

.221 

.699 

3 

2.374 

2.237 

.138 

.942 

5 

2.317 

2.239 

.258 

.897 

7 

2.913 

2.550 

.433 

.855 

9 

)  .641 

.893 

.748 

.344 

11 

1.800 

1.164 

637 

.647 

13 

2.o83 

1.877 

.806 

.700 

i: 

3.890 

3.331 

.559 

.856 

17 

2.760 

2.637 

.123 

.955 

19 

2.848 

1.540 

2.307 

.541 

21 

2.941 

1.499 

2.442 

.510 

2.760 

2.217 

.543 

.803 

25 

2.808 

1 . 584 

1.224 

.564 

The  dependency  of  8ext«  ®sca»  ^aba  ant*  M(0)  on  the  temperature 
is  not  expected  to  be  strong,  since  nx  and  n2  do  not  show  large 
variations  due  to  the  temperature.  As  will  be  seen  for  carbon 
particles,  relatively  large  variations  in  n2  due  to  temperature 
produce  only  moderate  variations  in  the  coefficients  and  negligible 
ones  in  the  scattering  function.  Similar  results  for  the  efficiency 
factors  of  single  aluminum  particles  were  also  found  by  Plass  (1965) . 


Figure  4-20.  Angular  Scattering  Phase  Function  M(6) 
for  A£203  at  X  =  1  pm. 


Figure  4-21.  Angular  Scattering  Phase  Function  M(e) 
for  A£203  at  X  ~  3  pm. 


Figure  4-22.  Angular  Scattering  Phase  Function  M(8) 
for  A£_0.  at  X  *  5  urn. 


Figure  4-23.  Angular  Scattering  Phase  Function  M(0) 
for  Ai203  at  X  =  7  yin. 


Figure  4-24.  Angular  Scattering  Phase  Function  M(0) 
for  A£20j  at  A  »  9  ym. 
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Figure  4-25.  Angular  Scattering  Phase  Function  M(6) 
for  A£-0.  at  X  =*  11  ym. 


Figure  4-29.  Angular  Scattering  Phase  Function  M(0) 
for  AijOj  at  X  <=  19  pin. 


Figure  4-30.  Angular  Scattering  Phase  Function  M(6) 
for  AfcjOj  at  X  *  21  pm. 


Figure  4-31.  Angular  Scattering  Phase  Function  M(6) 
for  Afc.O.  at  X  »  23  urn. 


Figure  4-32.  Angular  Scattering  Phase  Function  M(0) 
for  A£j03  at  X  =  25  pm. 


X  (put) 


Figure  4-33.  Ratio  of  Side  to  Forward  Scattering 
Components  vs.  Wavelength  for  A£«0~ 


For  the  second  particulate  species,  the  volume  coefficient 
and  angular  scattering  functions  have  been  calculated  for  carbon. 
The  discontinuous  power  law  was  used  to  represent  the  size  dis¬ 
tribution  as  reported  in  the  literature  (Rochelle,  1967),  viz., 


N(r)  -  N  for  . 01  <  r  <  . 03  pm 


N(r) 


for  . 03  <  r  <  2  pm 


The  total  number  of  particles  was  arbitrarily  chosen  to  be 
1013cm~3  which  corresponds  to  about  20%  mole  fraction 

(rmean  *  The  inPuts  for  amin'  amax’  m  and  n’ 

are  listed  in  Table  4-4.  In  order  to  study  the  influence 


TABLE  4-4.  INPUT  PARAMETERS  FOR  CARBON  PARTICLES 


A 

anln 

amax 

»1 

n2 

U000K) 

n  2 

(2500K) 

i 

.013 

18.6 

1.5 

.  3 

.5 

3 

.0042 

6.25 

2.7 

.37 

1.1 

5 

.0025 

3.  77 

3.0 

.4 

1.5 

7 

.0018 

2.69 

3.5 

.4 

1.9 

9 

.0014 

2.09 

3.7 

.4 

1.7 

11 

.0011 

1.71 

4.1 

.32 

1.3 

13 

.0010 

1.45 

4.4 

.29 

1.1 

15 

.0009 

1.26 

4.7 

.23 

.9 

17 

.0007 

1.11 

4.8 

.  21 

.75 

19 

.0006 

.99 

4.8 

.18 

.6 

21 

.0006 

.90 

4.8 

.  15 

.5 

23 

.0005 

.82 

4.6 

.  13 

.*■ 

Vr 

25 

.0005 

.75 

4.8 

.  12 

.3 
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of  the  temperature,  n2  is  given  for  1000  and  2500K.  The  results 
of  e  .  8  .  8  .  and  u>  for  1000  and  2500  K  are  shown  in 

Tables  4-5  and  4-6.  Although  the  absorption  index  is  higher  by  a 
factor  of  about  3  for  the  higher  temperature  throughout  the  wave¬ 
length  region  from  1  to  25  yra,  the  parameters  £Jext,  8sca,  Babs 
and  as  are  not  greatly  affected.  The  influence  on  the  scattering 
phase  function  is  even  less,  as  can  be  seen  in  Figures  4-34  through 
4-46 . 

As  for  A£20,,  we  have  used  the  ratio  of  the  phase  function 
at  90°  and  0°  to  indicate  the  degree  of  isotropic  scattering  of 
C(s) .  The  ratio  Is  plotted  in  Figure  4-47  and  compared  with  that 
of  A#.z03 . 


TABLE  4-5.  Results  for  Bext,  Bsca,  Babs 
and  a)  as  a  Function  of  X  for 
C(s)  Particles  in  Rocket 
Plume  Exhaust  (T  =  1000  K, 

N  =  10^-^/cm^. 


X(un) 

Bext 

®»ca 

B»ba 

w 

1 

1921 

936 

987 

.686 

3 

2136 

1107 

1027 

.519 

5 

2096 

1075 

1022 

.513 

7 

2006 

1039 

968 

.518 

9 

1865 

962 

922 

.510 

11 

1698 

812 

886 

.678 

13 

1578 

800 

778 

.507 

15 

1607 

717 

690 

.510 

17 

1305 

661 

663 

.507 

19 

1308 

693 

616 

.530 

21 

1229 

766 

68h 

.607 

23 

956 

588 

365 

.616 

25 

819 

675 

366 

.580 
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TABU  4-6.  Re.ults  for  B,xt.  08<,a,  6ab, 
and  u  as  a  Function  of  X  for 


C<8) 

Particles 

in 

Rocket 

Plume  Exhaust 

(T  - 

2500  K 

N  » 

lO^/cm-*) 

1<W»0 

*«xt 

«.c. 

_!ib* 

w 

1 

1927 

936 

992 

.486 

3 

2137 

1106 

1031 

.318 

s 

2093 

1082 

1013 

.524 

7 

1977 

1044 

932 

.528 

9 

1844 

927 

917 

.503 

U 

1698 

812 

886 

.478 

13 

1360 

714 

846 

.438 

13 

1413 

613 

803 

.431 

17 

1283 

320 

763 

.403 

19 

1198 

481 

717 

.402 

21 

1090 

461 

629 

.423 

23 

913 

397 

317 

.435 

23 

777 

327 

450 

.421 

i  •*»**« 


Figure  4-35.  Angular  Scattering  Phase  Function  M(0) 
for  Carbon  Particles  at  1  **  3  ytn. 
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Figure  4-36.  Angular  Scattering  Phase  Function  M(e) 
for  Carbon  Particles  at  X  =  5  (jm. 


Figure  4-37.  Angular  Scattering  Phase  Function  M(0) 
for  Carbon  Particles  at  *  -  7  um. 
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Figure  4-38*  Angular  Scattering  Phase  Function  M(0) 
for  Carbon  Particles  at  X  »  9  um. 
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Figure  4-39.  Angular  Scattering  Phase  Function  M(o) 
for  Carbon  Particles  at  A  =  11  urn. 
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Figure  4-41.  Angular  Scattering  Phase  Function  M(G) 
for  Carbon  Particles  at  A  -  15  »itn . 


Figure  4-42.  Angular  Scattering  Phase  Function  M(o) 
for  Carbon  Particles  at  A  =  17  iim. 
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Figure  4-43.  Angular  Scattering  Phase  Function  M(o) 
for  Carbon  Particles  at  A  =  19  ym. 
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-44.  Angular  Scattering  Phase  Function  M(e) 
for  Carbon  Particles  at  X  =  21pm. 


-45.  Angular  Scattering  Phase  Function  M(e) 
for  Carbon  Particles  at  X  *  23  pm. 
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-46.  Angular  Scattering  Phase  Function  M(o) 
for  Carbon  Particles  at  X  =25  pm. 
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Figure  4-47.  The  Ratio  of  the  Phase  Function  at 
90°  and  0°  Scattering  Angle  as  a 
Function  of  Aa  for  A£203  (+  —  +) 
and  C (s)  (o - o)  . 

4 . 5  Recommended  Data  File  and  Update  Procedure 

Based  on  the  studies  of  the  scattering  properties  of 
particles  performed  during  Subtask  I.A  it  is  recommended  that 
SIRRM  shall  include 


(1)  a  permanent  Data  File  of  ns(A,T),  ns(A,T) 

6ext(A’T)*  Sext(A>T)  and1Ms'd(A,o'T)  for 
each  species  s  and  for  each  size  distribution 
d;  (this  Data  File  may  be  minimized  to  as  few 
temperatures  and  size  distributions  as  indicated 
by  experience) ; 
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(2)  a  subroutine  that  allows  the  user  to  update  this 
Data  File  or  add  to  it  as  new  optical  para¬ 
meters  (n  (A,T),  n  O.T)  and  size  distributions) 

1  2 

become  available ; 

(3)  a  subroutine  that  compares  the  size  distribution 
as  calculated  by  the  JANNAF  flow  field  with  those 
stored  in  the  Data  File  and  that  will  indicate  to 

the  user  whether  or  not  the  stored  size  distribution(s) 
is  (are)  sufficient.  We  consider  also  the  option 
of  an  automatic  updating. 

The  advantages  of  the  recommended  approach  are  the  following: 

-  Computer  time  is  minimized  by  using  a  stored  Data 
File; 

•  Storage  requirements  are  minimized  by  using  as  few 
size  distributions  for  each  species  as  possible; 

•  Flexibility  is  maintained  by  allowing  the  user  to  up¬ 
date  the  Data  File  at  his  discretion; 

•  Indication  is  given  to  the  user  before  execution  whether 
an  update  of  the  Data  File  in  terms  of  new  and/or 
additional  size  distributions  is  necessary. 

A  tentative  logic  flow  diagram  of  the  recommended  approach  is 
shown  in  Figure  4-48.  After  the  JANNAF  flow  field  is  read  in  and 
the  relevant  flow  field  positions  for  the  line-of-sight  cal¬ 
culations  are  determined,  a  subroutine  calculates  the  analytical 
expression  for  the  size  distributions  for  each  position.  This 
is  necessary  since  the  JANNAF  flow  field  does  not  give  the 
analytical  expression  but  does  provide  the  number  density  for 
several  specific  particle  sizes  (private  communications,  S.  Dash, 

May  1979).  Different  analytical  expressions  will  be  used, 
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such  as  the  discontinuous  power  law  for  particles  with  very 
small  radii  (carhon)  and  modified  gamma  function  Cor  particles 

wllh  larger  rrnlil  (Av  0  ,  etc.) 

/  :•  \ 


The  computed  size  distributions  will  be  compared  with 
those  already  stored  in  the  Data  File.  If  the  differences 
are  minor  (criteria  must  yet  be  established)  no  update  is 
required.  If  the  differences  are  major,  3  ,  3__„  and  M(0) 

may  have  to  be  calculated  and  stored  in  the  Data  File.  If  the 
decision  for  an  update  is  desired  to  be  automatic  (probably  a 
default  input) ,  SIRRM  will  read  the  appropriate  n  and  n  and 
will  calculate  8  ,  3„„a  and  M(0).  If  the  decision  is  not  to 

be  automatic.,  the  user  has  the  option  whether  or  not  new  values 
of  8  ,  8  and  M(0)  should  be  calculated.  As  more  experience 

is  gained  about  the  influence  of  different  size  distributions 
on  the  radiation  field  using  the  six-flux  and  N-flux  models,  the 
decision  about  updating  the  Data  File  may  become  straightforward 
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